What can I do if I can’t install the third-party libraries from distribution packages?

What can I do if I can't install the third-party libraries from distribution packages?

Calculation of vertical velocity and preparation of output files

Calculation of the vertical velocity from the horizontal wind using the continuity equation

Calculation of the vertical velocity from the pre-calculated MARS parameter 77

Overview of the call of the``submit.py`` Python script and raw sequence of work steps in flex_extract.

what can i do if i can't install the third-party libraries from distribution packages?

installation problems with ``gateway`` and ``destination`` parameters

Fig. 1: Example of disaggregation scheme as implemented in older versions for an isolated precipitation event lasting one time interval (thick blue line). The amount of original precipitation after de-accumulation is given by the blue-shaded area. The green circles represent the discrete grid points after disaggregation and linearly interpolate in between them as indicated by the green line and the green-shaded area. Note that supporting points for the interpolation are shifted by half a time interval compared to the times when other meteorological fields are available (Hittmeir et al. 2018).

Fig. 2: Precipitation rate linearly interpolated using a sub-grid with two additional points. Colours as in Fig. 1 (Hittmeir et al. 2018).

Fig. 3: Schematic overview of the basic notation in a precipitation interval with the original precipitation rate g (green) as a step function and the interpolated data f (dark blue) as the piecewise linear function. The original time interval with fixed grid length \delta t is split equidistantly in three subintervals denoted by I_i^{1,2,3}, with the slopes in the subintervals as denoted by k_i^{1,2,3} . The sub-grid function values f_i, f_i^{1,2}, f_{i+1} are marked by red diamonds (Hittmeir et al. 2018).

calculation of vertical velocity and preparation of output files

calculation of the vertical velocity from the horizontal wind using the continuity equation

calculation of the vertical velocity from the pre-calculated mars parameter 77

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#*******************************************************************************
# @Author: Leopold Haimberger (University of Vienna)
#
# @Date: November 2015
#
# @Change History:
#
#    February 2018 - Anne Philipp (University of Vienna):
#        - applied PEP8 style guide
#        - added documentation
#        - moved install_args_and_control in here
#        - splitted code in smaller functions
#        - delete convert build files in here instead of compile job script
#        - changed static path names to variables from config file
#        - splitted install function into several smaller pieces
#        - use of tarfile package in python
#
# @License:
#    (C) Copyright 2014-2019.
#    Anne Philipp, Leopold Haimberger
#
#    SPDX-License-Identifier: CC-BY-4.0
#
#    This work is licensed under the Creative Commons Attribution 4.0
#    International License. To view a copy of this license, visit
#    http://creativecommons.org/licenses/by/4.0/ or send a letter to
#    Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
#
# @Methods:
#    main
#    get_install_cmdline_args
#    install_via_gateway
#    check_install_conditions
#    mk_tarball
#    un_tarball
#    mk_env_vars
#    mk_compilejob
#    mk_job_template
#    del_convert_build
#    mk_convert_build
#
'''This script installs the flex_extract program.

Depending on the selected installation environment (locally or on the
ECMWF server ecgate or cca) the program extracts the command line
arguments and the CONTROL file parameter and prepares the corresponding
environment.
The necessary files are collected in a tar ball and placed
at the target location. There, is is untared, the environment variables are
set, and the Fortran code is compiled.
If the ECMWF environment is selected, a job script is prepared and submitted
for the remaining configurations after putting the tar ball on the
target ECMWF server.

Type: install.py --help
to get information about command line parameters.
Read the documentation for usage instructions.
'''

# ------------------------------------------------------------------------------
# MODULES
from __future__ import print_function

import os
import sys
import subprocess
import tarfile
from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter

# software specific classes and modules from flex_extract
import _config
from Classes.ControlFile import ControlFile
from Classes.UioFiles import UioFiles
from Mods.tools import (make_dir, put_file_to_ecserver, submit_job_to_ecserver,
                        silent_remove, execute_subprocess, none_or_str)

# FUNCTIONS
def main():
    '''Controls the installation process. Calls the installation function
    if target is specified.

    Parameters
    ----------

    Return
    ------
    '''

    args = get_install_cmdline_args()
    c = ControlFile(args.controlfile)
    c.assign_args_to_control(args)
    check_install_conditions(c)

    if c.install_target.lower() != 'local': # ecgate or cca
        install_via_gateway(c)
    else: # local
        install_local(c)

    return

def get_install_cmdline_args():
    '''Decomposes the command line arguments and assigns them to variables.
    Apply default values for arguments not present.

    Return
    ------
    args : Namespace
        Contains the commandline arguments from script/program call.
    '''
    parser = ArgumentParser(description='Install flex_extract software '
                                        'locally or on ECMWF machines',
                            formatter_class=ArgumentDefaultsHelpFormatter)

    parser.add_argument('--target', dest='install_target',
                        type=none_or_str, default=None,
                        help="Valid targets: local | ecgate | cca , \
                        the latter two are at ECMWF")
    parser.add_argument("--makefile", dest="makefile",
                        help='Name of makefile for compiling the '
                        'Fortran program')
    parser.add_argument("--ecuid", dest="ecuid",
                        help='User id at ECMWF')
    parser.add_argument("--ecgid", dest="ecgid",
                        help='Group id at ECMWF')
    parser.add_argument("--gateway", dest="gateway",
                        help='Name of the local gateway server')
    parser.add_argument("--destination", dest="destination",
                        help='ecaccess association, e.g. '
                        'myUser@genericSftp')

    parser.add_argument("--installdir", dest="installdir",
                        help='Root directory of the '
                        'flex_extract installation')

    # arguments for job submission to ECMWF, only needed by submit.py
    parser.add_argument("--job_template", dest='job_template',
                        type=none_or_str, default="job.template",
                        help='Rudimentary template file to create a batch '
                        'job template for submission to ECMWF servers')

    parser.add_argument("--controlfile", dest="controlfile",
                        type=none_or_str, default='CONTROL_EA5',
                        help="A file that contains all CONTROL parameters.")

    args = parser.parse_args()

    return args

def install_via_gateway(c):
    '''Prepare data transfer to remote gateway and submit a job script which will
    install everything on the remote gateway.

    c : ControlFile
        Contains all the parameters of CONTROL file and
        command line.

    Return
    ------

    '''

    tarball_name = _config.FLEXEXTRACT_DIRNAME + '.tar'
    tar_file = os.path.join(_config.PATH_FLEXEXTRACT_DIR, tarball_name)

    mk_compilejob(c.makefile, c.ecuid, c.ecgid, c.installdir)

    mk_job_template(c.ecuid, c.ecgid, c.installdir)

    mk_env_vars(c.ecuid, c.ecgid, c.gateway, c.destination)

    mk_tarball(tar_file, c.install_target)

    put_file_to_ecserver(_config.PATH_FLEXEXTRACT_DIR, tarball_name,
                         c.install_target, c.ecuid, c.ecgid)

    submit_job_to_ecserver(c.install_target,
                           os.path.join(_config.PATH_REL_JOBSCRIPTS,
                                        _config.FILE_INSTALL_COMPILEJOB))

    silent_remove(tar_file)

    print('Job compilation script has been submitted to ecgate for ' +
          'installation in ' + c.installdir +
          '/' + _config.FLEXEXTRACT_DIRNAME)
    print('You should get an email with subject "flexcompile" within ' +
          'the next few minutes!')

    return

def install_local(c):
    '''Perform the actual installation on a local machine.

    Return
    ------

    '''

    tar_file = os.path.join(_config.PATH_FLEXEXTRACT_DIR,
                            _config.FLEXEXTRACT_DIRNAME + '.tar')

    if c.installdir == _config.PATH_FLEXEXTRACT_DIR:
        print('WARNING: installdir has not been specified')
        print('flex_extract will be installed in here by compiling the ' +
              'Fortran source in ' + _config.PATH_FORTRAN_SRC)
        os.chdir(_config.PATH_FORTRAN_SRC)
    else: # creates the target working directory for flex_extract
        c.installdir = os.path.expandvars(os.path.expanduser(
            c.installdir))
        if os.path.abspath(_config.PATH_FLEXEXTRACT_DIR) != \
           os.path.abspath(c.installdir):
            mk_tarball(tar_file, c.install_target)
            make_dir(os.path.join(c.installdir,
                                  _config.FLEXEXTRACT_DIRNAME))
            os.chdir(os.path.join(c.installdir,
            un_tarball(tar_file)
                                  _config.FLEXEXTRACT_DIRNAME,
                                  _config.PATH_REL_FORTRAN_SRC))

    # Create Fortran executable - CONVERT2
    print('Install ' +  _config.FLEXEXTRACT_DIRNAME + ' software at ' +
          c.install_target + ' in directory ' +
          os.path.abspath(c.installdir) + '\n')

    del_convert_build('.')
    mk_convert_build('.', c.makefile)

    os.chdir(_config.PATH_FLEXEXTRACT_DIR)
    if os.path.isfile(tar_file):
        os.remove(tar_file)

    return

def check_install_conditions(c):
    '''Checks necessary attributes and conditions
    for the installation, e.g. whether they exist and contain values.
    Otherwise set default values.

    Return
    ------

    '''

    if c.install_target and \
       c.install_target not in _config.INSTALL_TARGETS:
        print('ERROR: unknown or missing installation target ')
        print('target: ', c.install_target)
        print('please specify correct installation target ' +
              str(_config.INSTALL_TARGETS))
        print('use -h or --help for help')
        sys.exit(1)

    if c.install_target and c.install_target != 'local':
        if not c.ecgid or not c.ecuid:
            print('Please enter your ECMWF user id and group id '
                  ' with command line options --ecuid --ecgid')
            print('Try "' + sys.argv[0].split('/')[-1] + \
                  ' -h" to print usage information')
            print('Please consult ecaccess documentation or ECMWF user '
                  'support for further details.\n')
            sys.exit(1)
        if not c.gateway or not c.destination:
            print('WARNING: Parameters GATEWAY and DESTINATION were '
                  'not properly set for working on ECMWF server. \n'
                  'There will be no transfer of output files to the '
                  'local gateway server possible!')
        if not c.installdir:
            c.installdir = '${HOME}'
            c.installdir = _config.PATH_FLEXEXTRACT_DIR

    return

def mk_tarball(tarball_path, target):
    '''Creates a tarball with all necessary files which need to be sent to the
    installation directory.
    It does not matter whether this is local or remote.
    Collects all Python files, the Fortran source and makefiles,
    the ECMWF_ENV file, the CONTROL files as well as the
    template files.

    tarball_path : str
        The complete path to the tar file which will contain all
        relevant data for flex_extract.

    target : str
        The queue where the job is submitted to.

    Return
    ------

    '''

    print('Create tarball ...')

    # change to FLEXEXTRACT directory so that the tar can contain
    # relative pathes to the files and directories
    ecd = _config.PATH_FLEXEXTRACT_DIR + '/'
    os.chdir(ecd)

    # get lists of the files to be added to the tar file
    if target == 'local':
        ecmwf_env_file = []
        runfile = [os.path.relpath(x, ecd)
                   for x in UioFiles(_config.PATH_REL_RUN_DIR,
                                     'run_local.sh').files]
    else:
        ecmwf_env_file = [_config.PATH_REL_ECMWF_ENV]
                                     'run.sh').files]

    pyfiles = [os.path.relpath(x, ecd)
               for x in UioFiles(_config.PATH_REL_PYTHON_SRC, '*py').files]
    pytestfiles = [os.path.relpath(x, ecd)
                   for x in UioFiles(_config.PATH_REL_PYTHONTEST_SRC, '*py').files]
    controlfiles = [os.path.relpath(x, ecd)
                    for x in UioFiles(_config.PATH_REL_CONTROLFILES,
                                      'CONTROL*').files]
    testfiles = [os.path.relpath(x, ecd)
                 for x in UioFiles(_config.PATH_REL_TEST, '*').files]
    tempfiles = [os.path.relpath(x, ecd)
                 for x in UioFiles(_config.PATH_REL_TEMPLATES, '*.temp').files]
    nlfiles = [os.path.relpath(x, ecd)
               for x in UioFiles(_config.PATH_REL_TEMPLATES, '*.nl').files]
    gribtable = [os.path.relpath(x, ecd)
                 for x in UioFiles(_config.PATH_REL_TEMPLATES, '*grib*').files]
    ffiles = [os.path.relpath(x, ecd)
              for x in UioFiles(_config.PATH_REL_FORTRAN_SRC, '*.f90').files]
    hfiles = [os.path.relpath(x, ecd)
              for x in UioFiles(_config.PATH_REL_FORTRAN_SRC, '*.h').files]
    makefiles = [os.path.relpath(x, ecd)
                 for x in UioFiles(_config.PATH_REL_FORTRAN_SRC, 'makefile*').files]
    jobdir = [_config.PATH_REL_JOBSCRIPTS]

    # concatenate single lists to one for a better looping
    filelist = pyfiles + pytestfiles + controlfiles + tempfiles + nlfiles + \
               ffiles + gribtable + hfiles + makefiles + ecmwf_env_file + \
               runfile + jobdir + testfiles +\
               ['CODE_OF_CONDUCT.md', 'LICENSE.md', 'README.md']

    # create installation tar-file
    exclude_files = [".ksh", ".tar"]
    try:
        with tarfile.open(tarball_path, "w:gz") as tar_handle:
            for filename in filelist:
                tar_handle.add(filename, recursive=False,
                               filter=lambda tarinfo: None
                               if os.path.splitext(tarinfo.name)[1]
                               in exclude_files
                               else tarinfo)
    except tarfile.TarError as e:
        print('... ERROR: ' + str(e))

        sys.exit('\n... error occured while trying to create the tar-file ' +
                 str(tarball_path))

    return

def un_tarball(tarball_path):
    '''Extracts the given tarball into current directory.

    Return
    ------

    '''

    print('Untar ...')

    try:
        with tarfile.open(tarball_path) as tar_handle:
            tar_handle.extractall()
        sys.exit('\n... error occured while trying to read tar-file ' +
    except OSError as e:
        print('... ERROR CODE: ' + str(e.errno))
        print('... ERROR MESSAGE:\n \t ' + str(e.strerror))

    return

def mk_env_vars(ecuid, ecgid, gateway, destination):
    '''Creates a file named ECMWF_ENV which contains the
    necessary environmental variables at ECMWF servers.
    It is based on the template ECMWF_ENV.template.

    ecuid : str
        The user id on ECMWF server.

    ecgid : str
        The group id on ECMWF server.

    gateway : str
        The gateway server the user is using.

    destination : str
        The remote destination which is used to transfer files
        from ECMWF server to local gateway server.

    Return
    ------

    '''
    from genshi.template.text import NewTextTemplate
    from genshi.template import  TemplateLoader
    from genshi.template.eval import UndefinedError

    try:
        loader = TemplateLoader(_config.PATH_TEMPLATES, auto_reload=False)
        ecmwfvars_template = loader.load(_config.TEMPFILE_USER_ENVVARS,
                                         cls=NewTextTemplate)

        stream = ecmwfvars_template.generate(user_name=ecuid,
                                             user_group=ecgid,
                                             gateway_name=gateway,
                                             destination_name=destination
                                            )
    except UndefinedError as e:
        print('... ERROR ' + str(e))

        sys.exit('\n... error occured while trying to generate template ' +
                 _config.PATH_ECMWF_ENV)

    try:
        with open(_config.PATH_ECMWF_ENV, 'w') as f:
            f.write(stream.render('text'))

        sys.exit('\n... error occured while trying to write ' +

    return

def mk_compilejob(makefile, ecuid, ecgid, fp_root):
    '''Modifies the original job template file so that it is specified
    for the user and the environment were it will be applied. Result
    is stored in a new file "job.temp" in the python directory.

    makefile : str
        Name of the makefile which should be used to compile the Fortran
        program.

    ecuid : str

    ecgid : str

    fp_root : str
       Path to the root directory of FLEXPART environment or flex_extract
       environment.

    Return
    ------

    '''
    from genshi.template.eval import  UndefinedError

    if fp_root == '../':
        fp_root = '$HOME'

    try:
        compile_template = loader.load(_config.TEMPFILE_INSTALL_COMPILEJOB,
                                       cls=NewTextTemplate)

        stream = compile_template.generate(
            usergroup=ecgid,
            username=ecuid,
            version_number=_config._VERSION_STR,
            fp_root_scripts=fp_root,
            makefile=makefile,
            fortran_program=_config.FORTRAN_EXECUTABLE
        )

                 _config.TEMPFILE_INSTALL_COMPILEJOB)

    try:
        compilejob = os.path.join(_config.PATH_JOBSCRIPTS,
                                  _config.FILE_INSTALL_COMPILEJOB)

        with open(compilejob, 'w') as f:

                 compilejob)

    return

def mk_job_template(ecuid, ecgid, fp_root):
    is stored in a new file.

    ecuid : str

    ecgid : str

    Return
    ------

    '''

    fp_root_path_to_python = os.path.join(fp_root,
                                          _config.FLEXEXTRACT_DIRNAME,
                                          _config.PATH_REL_PYTHON_SRC)
    if '$' in fp_root_path_to_python:
        ind = fp_root_path_to_python.index('$')
        fp_root_path_to_python = fp_root_path_to_python[0:ind] + '$' + \
                                 fp_root_path_to_python[ind:]

    try:
        compile_template = loader.load(_config.TEMPFILE_INSTALL_JOB,

            fp_root_path=fp_root_path_to_python,
        )

                 _config.TEMPFILE_INSTALL_JOB)

    try:
        tempjobfile = os.path.join(_config.PATH_TEMPLATES,
                                   _config.TEMPFILE_JOB)

        with open(tempjobfile, 'w') as f:

                 tempjobfile)

    return

def del_convert_build(src_path):
    '''Clean up the Fortran source directory and remove all
    build files (e.g. \*.o, \*.mod and CONVERT2)

    src_path : str
        Path to the fortran source directory.

    Return
    ------

    '''

    modfiles = UioFiles(src_path, '*.mod')
    objfiles = UioFiles(src_path, '*.o')
    exefile = UioFiles(src_path, _config.FORTRAN_EXECUTABLE)

    modfiles.delete_files()
    objfiles.delete_files()
    exefile.delete_files()

    return

def mk_convert_build(src_path, makefile):
    '''Compiles the Fortran code and generates the executable.

        The name of the makefile which should be used.

    Return
    ------

    '''

    try:
        print('Using makefile: ' + makefile)
        p = subprocess.Popen(['make', '-f',
                              os.path.join(src_path, makefile)],
                             stdin=subprocess.PIPE,
                             stdout=subprocess.PIPE,
                             stderr=subprocess.PIPE,
                             bufsize=1)
        pout, perr = p.communicate()
        print(pout.decode())
        if p.returncode != 0:
            print(perr.decode())
            print('Please edit ' + makefile +
                  ' or try another makefile in the src directory.')
            print('Most likely GRIB_API_INCLUDE_DIR, GRIB_API_LIB '
                  'and EMOSLIB must be adapted.')
            print('Available makefiles:')
            print(UioFiles(src_path, 'makefile*'))
            sys.exit('Compilation failed!')
    except ValueError as e:
        print('ERROR: make of Fortran code failed:')
        print(e)
    else:
        execute_subprocess(['ls', '-l', 
                            os.path.join(src_path, _config.FORTRAN_EXECUTABLE)],
                           error_msg='FORTRAN EXECUTABLE COULD NOT BE FOUND!')

    return

if __name__ == "__main__":
    main()

# @Author: Anne Fouilloux (University of Oslo)
#
# @Date: October 2014
#
#
#    November 2015 - Leopold Haimberger (University of Vienna):
#        - job submission on ecgate and cca
#        - job templates suitable for twice daily operational dissemination
#
#        - minor changes in programming style (for consistence)
#        - changed path names to variables from config file
#        - added option for writing mars requests to extra file
#          additionally, as option without submitting the mars jobs
#        - splitted submit function to use genshi templates for the
#          job script and avoid code duplication
#
#
#
'''This script allows the user to extract meteorological fields from the ECMWF.

It prepares the settings for retrieving the data from ECMWF servers and
process the resulting files to prepare them for the use with FLEXPART or
FLEXTRA.

If it is supposed to work locally then it works through the necessary
functions get_mars_data and prepare_flexpart. Otherwise it prepares
a job script (korn shell) which will do the necessary work on the
ECMWF server. This script will de submitted via the ecaccess command
ecaccess-job-submit.

This file can also be imported as a module which then contains the following
functions:

    * main - the main function of the script
    * submit - calls mk_jobscript depending on operation mode and submits its
    * mk_jobscript - creates the job script from a template

Type: submit.py --help
'''

# MODULES

import os
import sys
from datetime import datetime, timedelta

from Mods.tools import (setup_controldata, normal_exit,
                        submit_job_to_ecserver)
from Mods.get_mars_data import get_mars_data
from Mods.prepare_flexpart import prepare_flexpart
#from Classes.ControlFile import ControlFile

# METHODS

    '''Get the arguments from script call and from CONTROL file.
    Decides from the argument "queue" if the local version
    is done "queue=None" or the gateway version with "queue=ecgate"
    or "queue=cca".

    Return
    ------

    '''

    c, ppid, queue, job_template = setup_controldata()

    # on local side
    # starting from an ECMWF server this would also be the local side
    called_from_dir = os.getcwd()
    if queue is None:
        if c.inputdir[0] != '/':
            c.inputdir = os.path.join(called_from_dir, c.inputdir)
        if c.outputdir[0] != '/':
            c.outputdir = os.path.join(called_from_dir, c.outputdir)
        get_mars_data(c)
        if c.request == 0 or c.request == 2:
            prepare_flexpart(ppid, c)
            exit_message = 'FLEX_EXTRACT IS DONE!'
        else:
            exit_message = 'PRINTING MARS_REQUESTS DONE!'
    # send files to ECMWF server
    else:
        submit(job_template, c, queue)
        exit_message = 'FLEX_EXTRACT JOB SCRIPT IS SUBMITED!'

    normal_exit(exit_message)

    return

def submit(jtemplate, c, queue):
    '''Prepares the job script and submits it to the specified queue.

    jtemplate : str
        Job template file from sub-directory "_templates" for
        submission to ECMWF. It contains all necessary
        module and variable settings for the ECMWF environment as well as
        the job call and mail report instructions.
        Default is "job.temp".

    queue : str
        Name of queue for submission to ECMWF (e.g. ecgate or cca )

    Return
    ------

    '''

    if not c.oper:
    # --------- create on demand job script ------------------------------------
        if c.purefc:
            print('---- Pure forecast mode! ----')
        else:
            print('---- On-demand mode! ----')

        job_file = os.path.join(_config.PATH_JOBSCRIPTS,
                                jtemplate[:-5] + '.ksh')

        # divide time periode into specified number of job chunks
        # to have multiple job scripts
        if c.job_chunk:
            start = datetime.strptime(c.start_date, '%Y%m%d')
            end = datetime.strptime(c.end_date, '%Y%m%d')
            chunk = timedelta(days=c.job_chunk)
            oneday = timedelta(days=1)

            while start <= end:
                if (start + chunk) <= end:
                    c.end_date = (start + chunk - oneday).strftime("%Y%m%d")
                else:
                    c.end_date = end.strftime("%Y%m%d")

                clist = c.to_list()
                mk_jobscript(jtemplate, job_file, clist)

                job_id = submit_job_to_ecserver(queue, job_file)
                print('The job id is: ' + str(job_id.strip()))

                start = start + chunk
                c.start_date = start.strftime("%Y%m%d")
        # submit a single job script
        else:
            clist = c.to_list()

            mk_jobscript(jtemplate, job_file, clist)

            job_id = submit_job_to_ecserver(queue, job_file)
            print('The job id is: ' + str(job_id.strip()))

    else:
    # --------- create operational job script ----------------------------------
        print('---- Operational mode! ----')

                                jtemplate[:-5] + 'oper.ksh')

        c.start_date = '${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}'
        c.end_date = '${MSJ_YEAR}${MSJ_MONTH}${MSJ_DAY}'
        c.basetime = '${MSJ_BASETIME}'
        if c.maxstep > 24:
            c.time = '${MSJ_BASETIME} {MSJ_BASETIME}'

        clist = c.to_list()

        mk_jobscript(jtemplate, job_file, clist)

        job_id = submit_job_to_ecserver(queue, job_file)
        print('The job id is: ' + str(job_id.strip()))

    print('You should get an email per job with subject flex.hostname.pid')

    return

def mk_jobscript(jtemplate, job_file, clist):
    '''Creates the job script from template.

    job_file : str
        Path to the job script file.

    clist : list of str
        Contains all necessary parameters for ECMWF CONTROL file.

    Return
    ------

    '''

    # load template and insert control content as list
    try:
        control_template = loader.load(jtemplate,

        stream = control_template.generate(control_content=clist)

        sys.exit('\n... error occured while trying to generate jobscript')

    # create jobscript file
    try:
        with open(job_file, 'w') as f:

        sys.exit('\n... error occured while trying to write ' + job_file)

    return

    main()

#
#
#
#   February 2018 - Anne Philipp (University of Vienna):
#        - applied some minor modifications in programming style/structure
#        - changed name of class Control to ControlFile for more
#          self-explanation naming
#        - outsource of class ControlFile
#        - initialisation of class attributes ( to avoid high number of
#          conditional statements and set default values )
#        - divided assignment of attributes and the check of conditions
#        - outsourced the commandline argument assignments to control attributes
#
#
#

# MODULES

import os
import sys

#pylint: disable=wrong-import-position
sys.path.append('../')
from Mods.tools import my_error
from Mods.checks import (check_grid, check_area, check_levels, check_purefc,
                         check_step, check_mail, check_queue, check_pathes,
                         check_dates, check_maxstep, check_type, check_request,
                         check_basetime, check_public, check_acctype,
                         check_acctime, check_accmaxstep, check_time,
                         check_logicals_type, check_len_type_time_step,
                         check_addpar, check_job_chunk, check_number)
#pylint: enable=wrong-import-position

# CLASS
class ControlFile(object):
    '''
    Contains the information which are stored in the CONTROL files.

    The CONTROL file is the steering part of the FLEXPART extraction
    software. All necessary parameters needed to retrieve the data fields
    from the MARS archive for driving FLEXPART are set in a CONTROL file.
    Some specific parameters like the start and end dates can be overwritten
    by the command line parameters, but in generall all parameters needed
    for a complete set of fields for FLEXPART can be set in the CONTROL file.

    Attributes
    controlfile : str
        The name of the control file to be processed. Default value is the
        filename passed to the init function when initialised.

    start_date : str
        The first day of the retrieval period. Default value is None.

    end_date :str
        The last day of the retrieval period. Default value is None.

    date_chunk : int
        Length of period for a single mars retrieval. Default value is 3.

    dtime :str
        The time step in hours. Default value is None.

    basetime : int
        The time for a half day retrieval. The 12 hours upfront are to be
        retrieved. Default value is None.

    maxstep : int
        The maximum forecast step for non flux data. Default value is None.

    type : list of str
        List of field type per retrieving hour. Default value is None.

    time : list of str
        List of retrieving times in hours. Default valuer is None.

    step : list of str or str
        List of forecast time steps in hours for non flux data.
        Default value is None.

    acctype : str
        The field type for the accumulated forecast fields.

    acctime : str
        The starting time of the accumulated forecasts. Default value is None.

    accmaxstep : int
        The maximum forecast step for the accumulated forecast fields
        (flux data). Default value is None.

    marsclass : str
        Characterisation of dataset. Default value is None.

    dataset : str
        For public datasets there is the specific naming and parameter
        dataset which has to be used to characterize the type of
        data. Default value is None.

    stream : str
        Identifies the forecasting system used to generate the data.

    number : str
        Selects the member in ensemble forecast run. Default value is 'OFF'.

    expver : str
        The version number of the dataset. Default value is '1'.

    gaussian : str
        This parameter is deprecated and should no longer be used.
        Specifies the desired type of Gaussian grid for the output.
        Default value is an empty string ''.

    grid : str
        Specifies the output grid which can be either a Gaussian grid
        or a Latitude/Longitude grid. Default value is None.

    area : str
        Specifies the desired sub-area of data to be extracted.

    left : str
        The western most longitude of the area to be extracted.

    lower : str
        The southern most latitude of the area to be extracted.

    upper : str
        The northern most latitued of the area to be extracted.

    right : str
        The eastern most longitude of the area to be extracted.

    level : str
        Specifies the maximum level. Default value is None.

    levelist : str
        Specifies the required level list. Default value is None.

    resol : str
        Specifies the desired triangular truncation of retrieved data,
        before carrying out any other selected post-processing.

    gauss : int
        Switch to select gaussian fields (1) or regular lat/lon (0).
        Default value is 0.

    accuracy : int
        Specifies the number of bits per value to be used in the
        generated GRIB coded fields. Default value is 24.

    omega : int
       Switch to select omega retrieval (1) or not (0). Default value is 0.

    omegadiff : int
        Switch to decide to calculate Omega and Dps/Dt from continuity
        equation for diagnostic purposes (1) or not (0). Default value is 0.

    eta : int
        Switch to select direct retrieval of etadot from MARS (1) or
        wether it has to be calculated (0). Then Default value is 0.

    etadiff : int
        Switch to select calculation of etadot and Dps/Dt from continuity

    etapar : int
        GRIB parameter Id for etadot fields. Default value is 77.

    dpdeta : int
        Switch to select multiplication of etadot with dpdeta.
        Default value is 1.

    smooth : int
        Spectral truncation of ETADOT after calculation on Gaussian grid.

    format : str
        The format of the GRIB data. Default value is 'GRIB1'.

    addpar : str
        List of additional surface level ECMWF parameter to be retrieved.

    prefix : str
        Prefix string for the final FLEXPART/FLEXTRA ready input files.
        Default value is 'EN'.

    cwc : int
        Switch to select wether the sum of cloud liquid water content and
        cloud ice water content should be retrieved. Default value is 0.

    wrf : int
        Switch to select further parameters for retrievment to support
        WRF simulations. Default value is 0.

    ecfsdir : str
        Path to the ECMWF storage  'ectmp:/${USER}/econdemand/'

    mailfail : list of str
        Email list for sending error log files from ECMWF servers.
        The email addresses should be seperated by a comma.
        Default value is ['${USER}'].

    mailops : list of str
        Email list for sending operational log files from ECMWF servers.

    ecstorage : int
        Switch to select storage of FLEXPART ready output files
        in the ECFS file system. Default value is 0.

    ectrans : int
        Switch to select the transfer of FLEXPART ready output files
        to the gateway server. Default value is 0.

    inputdir : str
        Path to the temporary directory for the retrieval grib files and
        other processing files. Default value is _config.PATH_INPUT_DIR.

    outputdir : str
        Path to the final directory where the final FLEXPART ready input
        files are stored. Default value is None.

    flexextractdir : str
        Path to the flex_extract root directory. Default value is
        _config.PATH_FLEXEXTRACT_DIR.

    exedir : str
        Path to the FORTRAN executable file. Default value is
        _config.PATH_FORTRAN_SRC.

    installdir : str
        Path to a FLEXPART root directory. Default value is None.

        Name of the makefile to be used for the Fortran program.
        Default value is 'Makefile.gfortran'.

        from ECMWF server to local gateway server. Default value is None.

        The gateway server the user is using. Default value is None.

    ecuid : str
        The user id on ECMWF server. Default value is None.

    ecgid : str
        The group id on ECMWF server. Default value is None.

    install_target : str
        Defines the location where the installation is to be done.

    debug : int
        Switch to keep temporary files at the end of postprocessing (1) or
        to delete all temporary files except the final output files (0).

    oper : int
        Switch to prepare the operational job script. Start date, end date and
        basetime will be prepared with environment variables.

    request : int
        Switch to select between just retrieving the data (0), writing the mars
        parameter values to a csv file (1) or doing both (2).

    public : int
        Switch to select kind of ECMWF Web Api access and the
        possible data sets. Public data sets (1) and Memberstate data sets (0).

    ec_api : boolean
        Tells wether the ECMWF Web API was able to load or not.

    cds_api : boolean
        Tells wether the CDS API was able to load or not.

    purefc : int
        Switch to decide wether the job is a pure forecast retrieval or
        coupled with analysis data. Default value is 0.

    rrint : int
        Switch to select between old precipitation disaggregation method (0)
        or the new IA3 disaggegration method (1). Default value is 0.

    doubleelda : int
        Switch to select the calculation of extra ensemble members for the
        ELDA stream. It doubles the amount of retrieved ensemble members.

    logicals : list of str
        List of the names of logical switches which controls the flow
        of the program. Default list is ['gauss', 'omega', 'omegadiff', 'eta',
        'etadiff', 'dpdeta', 'cwc', 'wrf', 'ecstorage',
        'ectrans', 'debug', 'request', 'public', 'purefc', 'rrint', 'doubleelda']
    '''

    def __init__(self, filename):
        '''Initialises the instance of ControlFile class and defines
        all class attributes with default values. Afterwards calls
        function __read_controlfile__ to read parameter from Control file.

        Parameters
        ----------
        filename : str
            Name of CONTROL file.

        Return
        ------

        '''

        # list of all possible class attributes and their default values
        self.controlfile = filename
        self.start_date = None
        self.end_date = None
        self.date_chunk = 3
        self.job_chunk = None
        self.dtime = None
        self.basetime = None
        self.maxstep = None
        self.type = None
        self.time = None
        self.step = None
        self.acctype = None
        self.acctime = None
        self.accmaxstep = None
        self.marsclass = None
        self.dataset = None
        self.stream = None
        self.number = 'OFF'
        self.expver = '1'
        self.gaussian = ''
        self.grid = None
        self.area = ''
        self.left = None
        self.lower = None
        self.upper = None
        self.right = None
        self.level = None
        self.levelist = None
        self.resol = None
        self.gauss = 0
        self.accuracy = 24
        self.omega = 0
        self.omegadiff = 0
        self.eta = 0
        self.etadiff = 0
        self.etapar = 77
        self.dpdeta = 1
        self.smooth = 0
        self.format = 'GRIB1'
        self.addpar = None
        self.prefix = 'EN'
        self.cwc = 0
        self.wrf = 0
        self.ecfsdir = 'ectmp:/${USER}/econdemand/'
        self.mailfail = ['${USER}']
        self.mailops = ['${USER}']
        self.ecstorage = 0
        self.ectrans = 0
        self.inputdir = _config.PATH_INPUT_DIR
        self.outputdir = None
        self.flexextractdir = _config.PATH_FLEXEXTRACT_DIR
        self.exedir = _config.PATH_FORTRAN_SRC
        self.installdir = None
        self.makefile = 'Makefile.gfortran'
        self.destination = None
        self.gateway = None
        self.ecuid = None
        self.ecgid = None
        self.install_target = None
        self.debug = 0
        self.oper = 0
        self.request = 0
        self.public = 0
        self.ec_api = None
        self.cds_api = None
        self.purefc = 0
        self.rrint = 0
        self.doubleelda = 0

        self.logicals = ['gauss', 'omega', 'omegadiff', 'eta', 'etadiff',
                         'dpdeta', 'cwc', 'wrf', 'ecstorage',
                         'ectrans', 'debug', 'oper', 'request', 'public',
                         'purefc', 'rrint', 'doubleelda']

        self._read_controlfile()

        return

    def _read_controlfile(self):
        '''Read CONTROL file and assign all CONTROL file variables.

        Return
        ------

        '''

        try:
            cfile = os.path.join(_config.PATH_CONTROLFILES, self.controlfile)
            with open(cfile) as f:
                fdata = f.read().split('\n')
        except IOError:
            print('Could not read CONTROL file "' + cfile + '"')
            print('Either it does not exist or its syntax is wrong.')
                      ' -h" to print usage information')

        # go through every line and store parameter
        for ldata in fdata:
            if ldata and ldata[0] == '#':
                # ignore comment line in control file
                continue
            if '#' in ldata:
                # cut off comment
                ldata = ldata.split('#')[0]
            data = ldata.split()
            if len(data) > 1:
                if 'm_' in data[0].lower():
                    data[0] = data[0][2:]
                if data[0].lower() == 'class':
                    data[0] = 'marsclass'
                if data[0].lower() == 'day1':
                    data[0] = 'start_date'
                if data[0].lower() == 'day2':
                    data[0] = 'end_date'
                if len(data) == 2:
                    if '$' in data[1]:
                        setattr(self, data[0].lower(), data[1])
                        while '$' in data[1]:
                            i = data[1].index('$')
                            j = data[1].find('{')
                            k = data[1].find('}')
                            var = os.getenv(data[1][j+1:k])
                            if var is not None:
                                data[1] = data[1][:i] + var + data[1][k+1:]
                            else:
                                my_error('Could not find variable '
                                         + data[1][j+1:k] + ' while reading ' +
                                         self.controlfile)
                        setattr(self, data[0].lower() + '_expanded', data[1])
                    else:
                        if data[1].lower() != 'none':
                            setattr(self, data[0].lower(), data[1])
                        else:
                            setattr(self, data[0].lower(), None)
                elif len(data) > 2:
                    setattr(self, data[0].lower(), (data[1:]))
            else:
                pass

        return

    def __str__(self):
        '''Prepares a string which have all the ControlFile class attributes
        with its associated values. Each attribute is printed in one line and
        in alphabetical order.

        Example
        -------
        'age': 10
        'color': 'Spotted'
        'kids': 0
        'legs': 2
        'name': 'Dog'
        'smell': 'Alot'

        Return
        ------
        string
            Single string of concatenated ControlFile class attributes
            with their values
        '''
        import collections

        attrs = vars(self).copy()
        attrs = collections.OrderedDict(sorted(attrs.items()))

        return '\n'.join("%s: %s" % item for item in attrs.items())

    def assign_args_to_control(self, args):
        '''Overwrites the existing ControlFile instance attributes with
        the command line arguments.

        args : Namespace
            Contains the commandline arguments from script/program call.

        Return
        ------

        '''

        # get dictionary of command line parameters and eliminate all
        # parameters which are None (were not specified)
        args_dict = vars(args)
        arguments = {k : args_dict[k] for k in args_dict
                     if args_dict[k] != None}

        # assign all passed command line arguments to ControlFile instance
        for k, v in arguments.items():
            setattr(self, str(k), v)

        return

    def assign_envs_to_control(self, envs):
        '''Assigns the ECMWF environment parameter.

        envs : dict of str
            Contains the ECMWF environment parameternames "ECUID", "ECGID",
            "DESTINATION" and "GATEWAY" with its corresponding values.
            They were read from the file "ECMWF_ENV".

        Return
        ------

        '''

        for k, v in envs.items():
            setattr(self, str(k).lower(), str(v))

        return

    def check_conditions(self, queue):
        '''Checks a couple of necessary attributes and conditions,
        such as if they exist and contain values.
        Otherwise set default values.

        queue : str
            Name of the queue if submitted to the ECMWF servers.
            Used to check if ecuid, ecgid, gateway and destination
            are set correctly and are not empty.

        Return
        ------

        '''
        check_logicals_type(self, self.logicals)

        self.mailfail = check_mail(self.mailfail)

        self.mailops = check_mail(self.mailops)

        check_queue(queue, self.gateway, self.destination,
                    self.ecuid, self.ecgid)

        self.outputdir, self.installdir = check_pathes(self.inputdir,
                                                       self.outputdir,
                                                       self.installdir,
                                                       self.flexextractdir)

        self.start_date, self.end_date = check_dates(self.start_date,
                                                     self.end_date)

        self.basetime = check_basetime(self.basetime)

        self.levelist, self.level = check_levels(self.levelist, self.level)

        self.step = check_step(self.step)

        self.maxstep = check_maxstep(self.maxstep, self.step)

        check_request(self.request,
                      os.path.join(self.inputdir, _config.FILE_MARS_REQUESTS))

        check_public(self.public, self.dataset)

        self.type = check_type(self.type, self.step)

        self.time = check_time(self.time)

        self.purefc = check_purefc(self.type)

        self.type, self.time, self.step = check_len_type_time_step(self.type,
                                                                   self.time,
                                                                   self.step,
                                                                   self.maxstep,
                                                                   self.purefc)

        self.acctype = check_acctype(self.acctype, self.type)

        self.acctime = check_acctime(self.acctime, self.marsclass,
                                     self.purefc, self.time)

        self.accmaxstep = check_accmaxstep(self.accmaxstep, self.marsclass,
                                           self.purefc, self.maxstep)

        self.grid = check_grid(self.grid)

        self.area = check_area(self.grid, self.area, self.upper, self.lower,
                               self.left, self.right)

        self.addpar = check_addpar(self.addpar)

        self.job_chunk = check_job_chunk(self.job_chunk)

        self.number = check_number(self.number)

        return

    def to_list(self):
        '''Just generates a list of strings containing the attributes and
        assigned values except the attributes "_expanded", "exedir",
        "flexextractdir" and "installdir".

        Return
        ------
        l : list of *
            A sorted list of the all ControlFile class attributes with
            their values except the attributes "_expanded", "exedir",
            "flexextractdir" and "installdir".
        '''

        attrs = collections.OrderedDict(sorted(vars(self).copy().items()))

        l = list()

        for item in attrs.items():
            if '_expanded' in item[0]:
                pass
            elif 'exedir' in item[0]:
                pass
            elif 'installdir' in item[0]:
                pass
            elif 'flexextractdir' in item[0]:
                pass
            else:
                if isinstance(item[1], list):
                    stot = ''
                    for s in item[1]:
                        stot += s + ' '

                    l.append("%s %s\n" % (item[0], stot))
                else:
                    l.append("%s %s\n" % item)

        return sorted(l)

#
#
#
#        - extended with class Control
#        - removed functions mkdir_p, daterange, years_between, months_between
#        - added functions darain, dapoly, to_param_id, init128, normal_exit,
#          my_error, clean_up, install_args_and_control,
#          interpret_args_and_control,
#        - removed function __del__ in class EIFLexpart
#        - added the following functions in EIFlexpart:
#            - create_namelist
#            - process_output
#            - deacc_fluxes
#        - modified existing EIFlexpart - functions for the use in
#          flex_extract
#        - retrieve also longer term forecasts, not only analyses and
#          short term forecast data
#        - added conversion into GRIB2
#
#        - removed function getFlexpartTime in class EcFlexpart
#        - outsourced class ControlFile
#        - outsourced class MarsRetrieval
#        - changed class name from EIFlexpart to EcFlexpart
#        - applied minor code changes (style)
#        - removed "dead code" , e.g. retrieval of Q since it is not needed
#        - removed "times" parameter from retrieve-method since it is not used
#        - seperated function "retrieve" into smaller functions (less code
#          duplication, easier testing)
#
#
#
#pylint: disable=unsupported-assignment-operation
# this is disabled because for this specific case its an error in pylint
#pylint: disable=consider-using-enumerate
# this is not useful in this case
#pylint: disable=unsubscriptable-object
# this error is a bug
#pylint: disable=ungrouped-imports
# not necessary that we group the imports
# MODULES

import os
import sys
import glob
import shutil

from Classes.GribUtil import GribUtil
from Mods.tools import (init128, to_param_id, silent_remove, product,
                        my_error, get_informations, get_dimensions,
                        execute_subprocess, to_param_id_with_tablenumber,
                        generate_retrieval_period_boundary)
from Classes.MarsRetrieval import MarsRetrieval
import Mods.disaggregation as disaggregation
# CLASS
class EcFlexpart(object):
    '''
    Class to represent FLEXPART specific ECMWF data.

    FLEXPART needs grib files in a specifc format. All necessary data fields
    for one time step are stored in a single file. The class represents an
    instance with all the parameter and settings necessary for retrieving
    MARS data and modifing them so they are fitting FLEXPART needs. The class
    is able to disaggregate the fluxes and convert grid types to the one needed
    by FLEXPART, therefore using the FORTRAN program.

    mreq_count : int
        Counter for the number of generated mars requests.

        Path to the directory where the retrieved data is stored.

        data.

        retrieved.

    dtime : str
        Time step in hours.

        The starting time from the accumulated forecasts.

    accmaxstep : str
        The maximum forecast step for the accumulated forecast fields.

        Characterisation of dataset.

        Selects the member in ensemble forecast run.

    resol : str

    accuracy : str
        generated GRIB coded fields.

        List of additional parameters to be retrieved.

    level : str
        Specifies the maximum level.

        The version of the dataset.

        Specifies the required levels.

    glevelist : str
        Specifies the required levels for gaussian grids.

    grid : str
        or a Latitude/Longitude grid.

    area : str

        Switch for definition of pure forecast mode or not.

    outputfilelist : list of str
        The final list of FLEXPART ready input files.

    types : dictionary
        Determines the combination of type of fields, time and forecast step
        to be retrieved.

    params : dictionary
        Collection of grid types and their corresponding parameters,
        levels, level types and the grid definition.

    server : ECMWFService or ECMWFDataServer
        This is the connection to the ECMWF data servers.

        Decides which Web API Server version is used.

    dates : str
        Contains start and end date of the retrieval in the format
        "YYYYMMDD/to/YYYYMMDD"
    '''

    # --------------------------------------------------------------------------
    # CLASS FUNCTIONS
    def __init__(self, c, fluxes=False):
        '''Creates an object/instance of EcFlexpart with the associated
        settings of its attributes for the retrieval.

        Parameters:
        -----------
        c : ControlFile
            Contains all the parameters of CONTROL file and
            command line.

        fluxes : boolean, optional
            Decides if the flux parameter settings are stored or
            the rest of the parameter list.
            Default value is False.

        Return
        ------

        '''
        # set a counter for the number of generated mars requests
        self.mreq_count = 0

        self.inputdir = c.inputdir
        self.dataset = c.dataset
        self.basetime = c.basetime
        self.dtime = c.dtime
        self.acctype = c.acctype
        self.acctime = c.acctime
        self.accmaxstep = c.accmaxstep
        self.marsclass = c.marsclass
        self.stream = c.stream
        self.number = c.number
        self.resol = c.resol
        self.accuracy = c.accuracy
        self.addpar = c.addpar
        self.level = c.level
        self.expver = c.expver
        self.levelist = c.levelist
        self.glevelist = '1/to/' + c.level # in case of gaussian grid
        self.gaussian = c.gaussian
        self.grid = c.grid
        self.area = c.area
        self.purefc = c.purefc
        self.outputfilelist = []

        # Define the different types of field combinations (type, time, step)
        self.types = {}
        # Define the parameters and their level types, level list and
        # grid resolution for the retrieval job
        self.params = {}

        if fluxes:
            self._create_params_fluxes()
        else:
            self._create_params(c.gauss, c.eta, c.omega, c.cwc, c.wrf)

        if fluxes:# and not c.purefc:
            self._create_field_types_fluxes()
        else:
            self._create_field_types(c.type, c.time, c.step)
        return

    def _create_field_types(self, ftype, ftime, fstep):
        '''Create the combination of field type, time and forecast step.

        ftype : list of str
            List of field types.

        ftime : list of str
            The time in hours of the field.

        fstep : str
            Specifies the forecast time step from forecast base time.
            Valid values are hours (HH) from forecast base time.

        Return
        ------

        '''
        i = 0
        for ty, st, ti in zip(ftype, fstep, ftime):
            btlist = list(range(len(ftime)))
            if self.basetime == 12:
                btlist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
            if self.basetime == 0:
                btlist = [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0]

            # if ((ty.upper() == 'AN' and (int(c.time[i]) % int(c.dtime)) == 0) or
                # (ty.upper() != 'AN' and (int(c.step[i]) % int(c.dtime)) == 0 and
                 # (int(c.step[i]) % int(c.dtime) == 0)) ) and \
                 # (int(c.time[i]) in btlist or c.purefc):

            if (i in btlist) or self.purefc:

                if ((ty.upper() == 'AN' and (int(ti) % int(self.dtime)) == 0) or
                    (ty.upper() != 'AN' and (int(st) % int(self.dtime)) == 0)):

                    if ty not in self.types.keys():
                        self.types[ty] = {'times': '', 'steps': ''}

                    if ti not in self.types[ty]['times']:
                        if self.types[ty]['times']:
                            self.types[ty]['times'] += '/'
                        self.types[ty]['times'] += ti

                    if st not in self.types[ty]['steps']:
                        if self.types[ty]['steps']:
                            self.types[ty]['steps'] += '/'
                        self.types[ty]['steps'] += st
            i += 1

        return

    def _create_field_types_fluxes(self):
        '''Create the combination of field type, time and forecast step
        for the flux data.

        Return
        ------

        '''
        if self.purefc:
            # need to retrieve forecasts for step 000 in case of pure forecast
            steps = '{}/to/{}/by/{}'.format(0, self.accmaxstep, self.dtime)
        else:
            steps = '{}/to/{}/by/{}'.format(self.dtime,
                                            self.accmaxstep,
                                            self.dtime)

        self.types[str(self.acctype)] = {'times': str(self.acctime),
                                         'steps': steps}
        return

    def _create_params(self, gauss, eta, omega, cwc, wrf):
        '''Define the specific parameter settings for retrievment.

        The different parameters need specific grid types and level types
        for retrievement. We might get following combination of types
        (depending on selection and availability):
        (These are short cuts for the grib file names (leading sequence)
        SH__ML, OG__ML, GG__ML, SH__SL, OG__SL, GG__SL, OG_OROLSM_SL
        where:
            SH = Spherical Harmonics, GG = Gaussian Grid, OG = Output Grid,
            ML = Model Level, SL = Surface Level

        For each of this combination there is a list of parameter names,
        the level type, the level list and the grid resolution.

        There are different scenarios for data extraction from MARS:
        1) Retrieval of etadot
           eta=1, gauss=0, omega=0
        2) Calculation of etadot from divergence
           eta=0, gauss=1, omega=0
        3) Calculation of etadot from omega (for makes sense for debugging)
           eta=0, gauss=0, omega=1
        4) Retrieval and Calculation of etadot (only for debugging)
           eta=1, gauss=1, omega=0
        5) Download also specific model and surface level data for FLEXPART-WRF

        gauss : int
            Gaussian grid is retrieved.

        eta : int
            Etadot parameter will be directly retrieved.

        omega : int
            The omega paramterwill be retrieved.

        cwc : int
            The cloud liquid and ice water content will be retrieved.

        wrf : int
            Additional model level and surface level data will be retrieved for
            WRF/FLEXPART-WRF simulations.

        Return
        ------

        '''
        # SURFACE FIELDS
        #-----------------------------------------------------------------------
        self.params['SH__SL'] = ['LNSP', 'ML', '1', 'OFF']
        self.params['OG__SL'] = ['SD/MSL/TCC/10U/10V/2T/2D/Z/LSM', \
                                 'SFC', '1', self.grid]
        if self.addpar:
            self.params['OG__SL'][0] += self.addpar

        if self.marsclass.upper() == 'EA' or self.marsclass.upper() == 'EP':
            self.params['OG_OROLSM__SL'] = ["SDOR/CVL/CVH/FSR",
                                            'SFC', '1', self.grid]
        else:
            self.params['OG_OROLSM__SL'] = ["SDOR/CVL/CVH/SR", \

        # MODEL LEVEL FIELDS
        self.params['OG__ML'] = ['T/Q', 'ML', self.levelist, self.grid]

        if not gauss and eta:
            self.params['OG__ML'][0] += '/U/V/ETADOT'
        elif gauss and not eta:
            self.params['GG__SL'] = ['Q', 'ML', '1',
                                     '{}'.format((int(self.resol) + 1) // 2)]
            self.params['SH__ML'] = ['U/V/D', 'ML', self.glevelist, 'OFF']
        elif not gauss and not eta:
            self.params['OG__ML'][0] += '/U/V'
        else:  # GAUSS and ETA
            print('Warning: Collecting etadot and parameters for gaussian grid '
                  'is a very costly parameter combination, '
                  'use this combination only for debugging!')
            self.params['GG__ML'] = ['U/V/D/ETADOT', 'ML', self.glevelist,

        if omega:
            self.params['OG__ML'][0] += '/W'

        if cwc:
            self.params['OG__ML'][0] += '/CLWC/CIWC'

        # ADDITIONAL FIELDS FOR FLEXPART-WRF MODEL (IF QUESTIONED)
        # ----------------------------------------------------------------------
        if wrf:
            # @WRF
            # THIS IS NOT YET CORRECTLY IMPLEMENTED !!!
            #
            # UNDER CONSTRUCTION !!!
            #

            print('WRF VERSION IS UNDER CONSTRUCTION!') # dummy argument

            #self.params['OG__ML'][0] += '/Z/VO'
            #if '/D' not in self.params['OG__ML'][0]:
            #    self.params['OG__ML'][0] += '/D'

            #wrf_sfc = ['SP','SKT','SST','CI','STL1','STL2', 'STL3','STL4',
            #           'SWVL1','SWVL2','SWVL3','SWVL4']
            #for par in wrf_sfc:
            #    if par not in self.params['OG__SL'][0]:
            #        self.params['OG__SL'][0] += '/' + par

        return

    def _create_params_fluxes(self):
        '''Define the parameter setting for flux data.

        Flux data are accumulated fields in time and are stored on the
        surface level. The leading short cut name for the grib files is:
        "OG_acc_SL" with OG for Regular Output Grid, SL for Surface Level, and
        acc for Accumulated Grid.
        The params dictionary stores a list of parameter names, the level type,
        the level list and the grid resolution.

        The flux data are: LSP/CP/SSHF/EWSS/NSSS/SSR

        Return
        ------

        '''
        self.params['OG_acc_SL'] = ["LSP/CP/SSHF/EWSS/NSSS/SSR",
                                    'SFC', '1', self.grid]
        return

    def _mk_targetname(self, ftype, param, date):
        '''Creates the filename for the requested grib data to be stored in.
        This name is passed as the "target" parameter in the request.

        ftype : str
            Shortcut name of the type of the field. E.g. AN, FC, PF, ...

        param : str
            Shortcut of the grid type. E.g. SH__ML, SH__SL, GG__ML,
            GG__SL, OG__ML, OG__SL, OG_OROLSM_SL, OG_acc_SL

        date : str
            The date period of the grib data to be stored in this file.

        Return
        ------
        targetname : str
            The target filename for the grib data.
        '''
        targetname = (self.inputdir + '/' + ftype + param + '.' + date + '.' +
                      str(os.getppid()) + '.' + str(os.getpid()) + '.grb')

        return targetname

    def _start_retrievement(self, request, par_dict):
        '''Creates the Mars Retrieval and prints or submits the request
        depending on the status of the request variable.

        request : int
            Selects the mode of retrieval.
            0: Retrieves the data from ECMWF.
            1: Prints the mars requests to an output file.
            2: Retrieves the data and prints the mars request.

        par_dict : dictionary
            Contains all parameter which have to be set for creating the
            Mars Retrievals. The parameter are:
            marsclass, dataset, stream, type, levtype, levelist, resol,
            gaussian, accuracy, grid, target, area, date, time, number,
            step, expver, param

        Return
        ------

        '''
        # increase number of mars requests
        self.mreq_count += 1

        MR = MarsRetrieval(self.server,
                           self.public,
                           marsclass=par_dict['marsclass'],
                           dataset=par_dict['dataset'],
                           stream=par_dict['stream'],
                           type=par_dict['type'],
                           levtype=par_dict['levtype'],
                           levelist=par_dict['levelist'],
                           resol=par_dict['resol'],
                           gaussian=par_dict['gaussian'],
                           accuracy=par_dict['accuracy'],
                           grid=par_dict['grid'],
                           target=par_dict['target'],
                           area=par_dict['area'],
                           date=par_dict['date'],
                           time=par_dict['time'],
                           number=par_dict['number'],
                           step=par_dict['step'],
                           expver=par_dict['expver'],
                           param=par_dict['param'])

        if request == 0:
            MR.display_info()
            MR.data_retrieve()
        elif request == 1:
            MR.print_infodata_csv(self.inputdir, self.mreq_count)
        elif request == 2:
        else:
            print('Failure')

        return

    def _mk_index_values(self, inputdir, inputfiles, keys):
        '''Creates an index file for a set of grib parameter keys.
        The values from the index keys are returned in a list.

        keys : dictionary
            List of parameter names which serves as index.

        inputfiles : UioFiles
            Contains a list of files.

        Return
        ------
        iid : codes_index
            This is a grib specific index structure to access
            messages in a file.

        index_vals : list of list  of str
            Contains the values from the keys used for a distinct selection
            of grib messages in processing  the grib files.
            Content looks like e.g.:
            index_vals[0]: ('20171106', '20171107', '20171108') ; date
            index_vals[1]: ('0', '1200', '1800', '600') ; time
            index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange
        '''
        from eccodes import codes_index_get

        iid = None
        index_keys = keys

        indexfile = os.path.join(inputdir, _config.FILE_GRIB_INDEX)
        silent_remove(indexfile)
        grib = GribUtil(inputfiles.files)
        # creates new index file
        iid = grib.index(index_keys=index_keys, index_file=indexfile)

        # read the values of index keys
        index_vals = []
        for key in index_keys:
            key_vals = codes_index_get(iid, key)
            # have to sort the key values for correct order,
            # therefore convert to int first
            key_vals = [int(k) for k in key_vals]
            key_vals.sort()
            key_vals = [str(k) for k in key_vals]
            index_vals.append(key_vals)
            # index_vals looks for example like:
            # index_vals[0]: ('20171106', '20171107', '20171108') ; date
            # index_vals[1]: ('0', '1200') ; time
            # index_vals[2]: (3', '6', '9', '12') ; stepRange

        return iid, index_vals

    def retrieve(self, server, dates, public, request, inputdir='.'):
        '''Finalizing the retrieval information by setting final details
        depending on grid type.
        Prepares MARS retrievals per grid type and submits them.

        server : ECMWFService or ECMWFDataServer
            The connection to the ECMWF server. This is different
            for member state users which have full access and non
            member state users which have only access to the public
            data sets. The decision is made from command line argument
            "public"; for public access its True (ECMWFDataServer)
            for member state users its False (ECMWFService)

        dates : str
            Contains start and end date of the retrieval in the format
            "YYYYMMDD/to/YYYYMMDD"

        inputdir : str, optional
            Path to the directory where the retrieved data is about
            to be stored. The default is the current directory ('.').

        Return
        ------

        '''
        self.dates = dates
        self.server = server
        self.public = public
        self.inputdir = inputdir
        oro = False

        # define times with datetime module
        t12h = timedelta(hours=12)
        t24h = timedelta(hours=24)

        # dictionary which contains all parameter for the mars request,
        # entries with a "None" will change in different requests and will
        # therefore be set in each request seperately
        retr_param_dict = {'marsclass':self.marsclass,
                           'dataset':self.dataset,
                           'stream':None,
                           'type':None,
                           'levtype':None,
                           'levelist':None,
                           'resol':self.resol,
                           'gaussian':None,
                           'accuracy':self.accuracy,
                           'grid':None,
                           'target':None,
                           'area':None,
                           'date':None,
                           'time':None,
                           'number':self.number,
                           'step':None,
                           'expver':self.expver,
                           'param':None}

        for ftype in sorted(self.types):
            # ftype contains field types such as
            #     [AN, FC, PF, CV]
            for pk, pv in sorted(self.params.items()):
                # pk contains one of these keys of params
                #     [SH__ML, SH__SL, GG__ML, GG__SL, OG__ML, OG__SL,
                #      OG_OROLSM_SL, OG_acc_SL]
                # pv contains all of the items of the belonging key
                #     [param, levtype, levelist, grid]
                if isinstance(pv, str):
                    continue
                retr_param_dict['type'] = ftype
                retr_param_dict['time'] = self.types[ftype]['times']
                retr_param_dict['step'] = self.types[ftype]['steps']
                retr_param_dict['date'] = self.dates
                retr_param_dict['stream'] = self.stream
                retr_param_dict['target'] = \
                    self._mk_targetname(ftype,
                                        pk,
                                        retr_param_dict['date'].split('/')[0])
                table128 = init128(_config.PATH_GRIBTABLE)
                ids = to_param_id_with_tablenumber(pv[0], table128)
                retr_param_dict['param'] = ids
                retr_param_dict['levtype'] = pv[1]
                retr_param_dict['levelist'] = pv[2]
                retr_param_dict['grid'] = pv[3]
                retr_param_dict['area'] = self.area
                retr_param_dict['gaussian'] = self.gaussian

                if pk == 'OG_OROLSM__SL' and not oro:
                    oro = True
                    # in CERA20C (class EP) there is no stream "OPER"!
                    if self.marsclass.upper() != 'EP':
                        retr_param_dict['stream'] = 'OPER'
                    retr_param_dict['type'] = 'AN'
                    retr_param_dict['time'] = '00'
                    retr_param_dict['step'] = '000'
                    retr_param_dict['date'] = self.dates.split('/')[0]
                    retr_param_dict['target'] = self._mk_targetname('',
                                                                    pk,
                                                                    retr_param_dict['date'])
                elif pk == 'OG_OROLSM__SL' and oro:
                    continue
                if pk == 'GG__SL' and pv[0] == 'Q':
                    retr_param_dict['area'] = ""
                    retr_param_dict['gaussian'] = 'reduced'
                if ftype.upper() == 'FC' and \
                    'acc' not in retr_param_dict['target']:
                    if (int(retr_param_dict['time'][0]) +
                        int(retr_param_dict['step'][0])) > 23:
                        dates = retr_param_dict['date'].split('/')
                        sdate = datetime.strptime(dates[0], '%Y%m%d%H')
                        sdate = sdate - timedelta(days=1)
                        retr_param_dict['date'] = '/'.join(
                            [sdate.strftime("%Y%m%d")] +
                            retr_param_dict['date'][1:])

                        print('CHANGED FC start date to ' +
                              sdate.strftime("%Y%m%d") +
                              ' to accomodate TIME=' +
                              retr_param_dict['time'][0] +
                              ', STEP=' +
                              retr_param_dict['time'][0])

    # ------  on demand path  --------------------------------------------------
                if self.basetime is None:
                    # ******* start retrievement
                    self._start_retrievement(request, retr_param_dict)
    # ------  operational path  ------------------------------------------------
                else:
                    # check if mars job requests fields beyond basetime.
                    # if yes eliminate those fields since they may not
                    # be accessible with user's credentials

                    enddate = retr_param_dict['date'].split('/')[-1]
                    elimit = datetime.strptime(enddate + str(self.basetime),
                                               '%Y%m%d%H')

                    if self.basetime == 12:
                        # --------------  flux data ----------------------------
                        if 'acc' in pk:
                            startdate = retr_param_dict['date'].split('/')[0]
                            enddate = datetime.strftime(elimit - t24h, '%Y%m%d')
                            retr_param_dict['date'] = '/'.join([startdate,
                                                                'to',
                                                                enddate])

                            # ******* start retrievement
                            self._start_retrievement(request, retr_param_dict)

                            retr_param_dict['date'] = \
                                datetime.strftime(elimit - t12h, '%Y%m%d')
                            retr_param_dict['time'] = '00'
                            retr_param_dict['target'] = \
                                self._mk_targetname(ftype, pk,
                                                    retr_param_dict['date'])

                        # --------------  non flux data ------------------------
                        else:

                    elif self.basetime == 0:
#                        retr_param_dict['date'] = \
#                            datetime.strftime(elimit - t24h, '%Y%m%d')

                        timesave = ''.join(retr_param_dict['time'])

                        if all(['/' in retr_param_dict['time'],
                                pk != 'OG_OROLSM__SL',
                                'acc' not in pk]):
                            times = retr_param_dict['time'].split('/')
                            steps = retr_param_dict['step'].split('/')

                            while int(times[0]) + int(steps[0]) <= 12:
                                times = times[1:]
                                if len(times) > 1:
                                    retr_param_dict['time'] = '/'.join(times)
                                else:
                                    retr_param_dict['time'] = times[0]

                        if all([pk != 'OG_OROLSM__SL',
                                int(retr_param_dict['step'].split('/')[0]) == 0,
                                int(timesave.split('/')[0]) == 0]):

                                datetime.strftime(elimit, '%Y%m%d')
                            retr_param_dict['step'] = '000'

                        # ******* start retrievement
                        self._start_retrievement(request, retr_param_dict)
                    else:
                        raise ValueError('ERROR: Basetime has an invalid value '
                                         '-> {}'.format(str(self.basetime)))

        if request == 0 or request == 2:
            print('MARS retrieve done ... ')
            print('MARS request printed ...')

        return

    def write_namelist(self, c):
        '''Creates a namelist file in the temporary directory and writes
        the following values to it: maxl, maxb, mlevel,
        mlevelist, mnauf, metapar, rlo0, rlo1, rla0, rla1,
        momega, momegadiff, mgauss, msmooth, meta, metadiff, mdpdeta

                Name of the namelist file.

        Return
        ------

        '''

        from genshi.template.text import NewTextTemplate
        from genshi.template import  TemplateLoader
        from genshi.template.eval import  UndefinedError
        import numpy as np

        try:
            loader = TemplateLoader(_config.PATH_TEMPLATES, auto_reload=False)
            namelist_template = loader.load(_config.TEMPFILE_NAMELIST,
                                            cls=NewTextTemplate)

            self.inputdir = c.inputdir
            area = np.asarray(self.area.split('/')).astype(float)
            grid = np.asarray(self.grid.split('/')).astype(float)

            if area[1] > area[3]:
                area[1] -= 360
            maxl = int(round((area[3] - area[1]) / grid[1])) + 1
            maxb = int(round((area[0] - area[2]) / grid[0])) + 1

            stream = namelist_template.generate(
                maxl=str(maxl),
                maxb=str(maxb),
                mlevel=str(self.level),
                mlevelist=str(self.levelist),
                mnauf=str(self.resol),
                metapar='77',
                rlo0=str(area[1]),
                rlo1=str(area[3]),
                rla0=str(area[2]),
                rla1=str(area[0]),
                momega=str(c.omega),
                momegadiff=str(c.omegadiff),
                mgauss=str(c.gauss),
                msmooth=str(c.smooth),
                meta=str(c.eta),
                metadiff=str(c.etadiff),
                mdpdeta=str(c.dpdeta)
            )
        except UndefinedError as e:
            print('... ERROR ' + str(e))

            sys.exit('\n... error occured while trying to generate namelist ' +
                     _config.TEMPFILE_NAMELIST)
        except OSError as e:
            print('... ERROR CODE: ' + str(e.errno))
            print('... ERROR MESSAGE:\n \t ' + str(e.strerror))

            sys.exit('\n... error occured while trying to generate template ' +

        try:
            namelistfile = os.path.join(self.inputdir, _config.FILE_NAMELIST)

            with open(namelistfile, 'w') as f:
                f.write(stream.render('text'))

            sys.exit('\n... error occured while trying to write ' +
                     namelistfile)

        return

    def deacc_fluxes(self, inputfiles, c):
        '''De-accumulate and disaggregate flux data.

        Goes through all flux fields in ordered time and de-accumulate
        the fields. Afterwards the fields are disaggregated in time.
        Different versions of disaggregation is provided for rainfall
        data (darain, modified linear) and the surface fluxes and
        stress data (dapoly, cubic polynomial).

            Contains the list of files that contain flux data.

        Return
        ------

        '''
        from eccodes import (codes_index_select, codes_get,
                             codes_get_values, codes_set_values, codes_set,
                             codes_write, codes_release, codes_new_from_index,
                             codes_index_release)

        table128 = init128(_config.PATH_GRIBTABLE)
        # get ids from the flux parameter names
        pars = to_param_id(self.params['OG_acc_SL'][0], table128)

        iid = None
        index_vals = None

        # get the values of the keys which are used for distinct access
        # of grib messages via product and save the maximum number of
        # ensemble members if there is more than one
        if '/' in self.number:
            # more than one ensemble member is selected
            index_keys = ["number", "date", "time", "step"]
            # maximum ensemble number retrieved
            # + 1 for the control run (ensemble number 0)
            maxnum = int(self.number.split('/')[-1]) + 1
            # remember the index of the number values
            index_number = index_keys.index('number')
            # empty set to save ensemble numbers which were already processed
            ens_numbers = set()
            # index for the ensemble number
            inumb = 0
        else:
            index_keys = ["date", "time", "step"]
            # maximum ensemble number
            maxnum = None

        # get sorted lists of the index values
        # this is very important for disaggregating
        # the flux data in correct order
        iid, index_vals = self._mk_index_values(c.inputdir,
                                                inputfiles,
                                                index_keys)
        # index_vals looks like e.g.:
        # index_vals[0]: ('20171106', '20171107', '20171108') ; date
        # index_vals[1]: ('0', '600', '1200', '1800') ; time
        # index_vals[2]: ('0', '3', '6', '9', '12') ; stepRange

        if c.rrint:
            # set start and end timestamps for retrieval period
            if not c.purefc:
                start_date = datetime.strptime(c.start_date + '00', '%Y%m%d%H')
                end_date = datetime.strptime(c.end_date + '23', '%Y%m%d%H')
            else:
                sdate_str = c.start_date + '{:0>2}'.format(index_vals[1][0])
                start_date = datetime.strptime(sdate_str, '%Y%m%d%H')
                edate_str = c.end_date + '{:0>2}'.format(index_vals[1][-1])
                end_date = datetime.strptime(edate_str, '%Y%m%d%H')
                end_date = end_date + timedelta(hours=c.maxstep)

            # get necessary grid dimensions from grib files for storing the
            # precipitation fields
            info = get_informations(os.path.join(c.inputdir,
                                                 inputfiles.files[0]))
            dims = get_dimensions(info, c.purefc, c.dtime, index_vals,
                                  start_date, end_date)

            # create empty numpy arrays
            if not maxnum:
                lsp_np = np.zeros((dims[1] * dims[0], dims[2]), dtype=np.float64)
                cp_np = np.zeros((dims[1] * dims[0], dims[2]), dtype=np.float64)
            else:
                lsp_np = np.zeros((maxnum, dims[1] * dims[0], dims[2]), dtype=np.float64)
                cp_np = np.zeros((maxnum, dims[1] * dims[0], dims[2]), dtype=np.float64)

            # index counter for time line
            it_lsp = 0
            it_cp = 0

            # store the order of date and step
            date_list = []
            step_list = []

        # initialize dictionaries to store flux values per parameter
        orig_vals = {}
        deac_vals = {}
        for p in pars:
            orig_vals[p] = []
            deac_vals[p] = []

        # "product" genereates each possible combination between the
        # values of the index keys
        for prod in product(*index_vals):
            # e.g. prod = ('20170505', '0', '12')
            #             ( date     ,time, step)

            print('CURRENT PRODUCT: ', prod)

            # the whole process has to be done for each seperate ensemble member
            # therefore, for each new ensemble member we delete old flux values
            # and start collecting flux data from the beginning time step
            if maxnum and prod[index_number] not in ens_numbers:
                ens_numbers.add(prod[index_number])
                inumb = len(ens_numbers) - 1
                # re-initialize dictionaries to store flux values per parameter
                # for the next ensemble member
                it_lsp = 0
                it_cp = 0
                orig_vals = {}
                deac_vals = {}
                for p in pars:
                    orig_vals[p] = []
                    deac_vals[p] = []

            for i in range(len(index_keys)):
                codes_index_select(iid, index_keys[i], prod[i])

            # get first id from current product
            gid = codes_new_from_index(iid)

            # if there is no data for this specific time combination / product
            # skip the rest of the for loop and start with next timestep/product
            if not gid:
                continue

            # create correct timestamp from the three time informations
            cdate = str(codes_get(gid, 'date'))
            time = codes_get(gid, 'time') // 100  # integer
            step = codes_get(gid, 'step') # integer
            ctime = '{:0>2}'.format(time)

            t_date = datetime.strptime(cdate + ctime, '%Y%m%d%H')
            t_dt = t_date + timedelta(hours=step)
            t_m1dt = t_date + timedelta(hours=step-int(c.dtime))
            t_m2dt = t_date + timedelta(hours=step-2*int(c.dtime))
            if c.basetime is not None:
                t_enddate = datetime.strptime(c.end_date + str(c.basetime),
                                              '%Y%m%d%H')
            else:
                t_enddate = t_date + timedelta(2*int(c.dtime))

            # if necessary, add ensemble member number to filename suffix
            # otherwise, add empty string
            if maxnum:
                numbersuffix = '.N{:0>3}'.format(int(prod[index_number]))
            else:
                numbersuffix = ''

            if c.purefc:
                fnout = os.path.join(c.inputdir, 'flux' +
                                     t_date.strftime('%Y%m%d.%H') +
                                     '.{:0>3}'.format(step-2*int(c.dtime)) +
                                     numbersuffix)
                gnout = os.path.join(c.inputdir, 'flux' +
                                     '.{:0>3}'.format(step-int(c.dtime)) +
                hnout = os.path.join(c.inputdir, 'flux' +
                                     '.{:0>3}'.format(step) +
            else:
                                     t_m2dt.strftime('%Y%m%d%H') + numbersuffix)
                                     t_m1dt.strftime('%Y%m%d%H') + numbersuffix)
                                     t_dt.strftime('%Y%m%d%H') + numbersuffix)

            print("outputfile = " + fnout)
            f_handle = open(fnout, 'wb')
            h_handle = open(hnout, 'wb')
            g_handle = open(gnout, 'wb')

            # read message for message and store relevant data fields, where
            # data keywords are stored in pars
            while True:
                if not gid:
                    break
                parId = codes_get(gid, 'paramId') # integer
                step = codes_get(gid, 'step') # integer
                time = codes_get(gid, 'time') # integer
                ni = codes_get(gid, 'Ni') # integer
                nj = codes_get(gid, 'Nj') # integer
                if parId not in orig_vals.keys():
                    # parameter is not a flux, find next one
                    continue

                # define conversion factor
                if parId == 142 or parId == 143:
                    fak = 1. / 1000.
                else:
                    fak = 3600.

                # get parameter values and reshape
                values = codes_get_values(gid)
                values = (np.reshape(values, (nj, ni))).flatten() / fak

                # save the original and accumulated values
                orig_vals[parId].append(values[:])

                if c.marsclass.upper() == 'EA' or step <= int(c.dtime):
                    # no de-accumulation needed
                    deac_vals[parId].append(values[:] / int(c.dtime))
                else:
                    # do de-accumulation
                    deac_vals[parId].append(
                        (orig_vals[parId][-1] - orig_vals[parId][-2]) /
                        int(c.dtime))

                # store precipitation if new disaggregation method is selected
                # only the exact days are needed
                if c.rrint:
                    if start_date <= t_dt <= end_date:
                        if not c.purefc:
                            if t_dt not in date_list:
                                date_list.append(t_dt)
                                step_list = [0]
                        else:
                            if t_date not in date_list:
                                date_list.append(t_date)
                            if step not in step_list:
                                step_list.append(step)
                        # store precipitation values
                        if maxnum and parId == 142:
                            lsp_np[inumb, :, it_lsp] = deac_vals[parId][-1][:]
                            it_lsp += 1
                        elif not maxnum and parId == 142:
                            lsp_np[:, it_lsp] = deac_vals[parId][-1][:]
                            it_lsp += 1
                        elif maxnum and parId == 143:
                            cp_np[inumb, :, it_cp] = deac_vals[parId][-1][:]
                            it_cp += 1
                        elif not maxnum and parId == 143:
                            cp_np[:, it_cp] = deac_vals[parId][-1][:]
                            it_cp += 1

                # information printout
                print(parId, time, step, len(values), values[0], np.std(values))

                # length of deac_vals[parId] corresponds to the
                # number of time steps, max. 4 are needed for disaggegration
                # with the old and original method
                # run over all grib messages and perform
                # shifting in time
                if len(deac_vals[parId]) >= 3:
                    if len(deac_vals[parId]) > 3:
                        if not c.rrint and (parId == 142 or parId == 143):
                            values = disaggregation.darain(deac_vals[parId])
                        else:
                            values = disaggregation.dapoly(deac_vals[parId])

                        if not (step == c.maxstep and c.purefc \
                                or t_dt == t_enddate):
                            # remove first time step in list to shift
                            # time line
                            orig_vals[parId].pop(0)
                            deac_vals[parId].pop(0)
                    else:
                        # if the third time step is read (per parId),
                        # write out the first one as a boundary value
                        if c.purefc:
                            values = deac_vals[parId][1]
                        else:
                            values = deac_vals[parId][0]

                    if not (c.rrint and (parId == 142 or parId == 143)):
                        codes_set_values(gid, values)

                            codes_set(gid, 'stepRange', max(0, step-2*int(c.dtime)))
                        else:
                            codes_set(gid, 'stepRange', 0)
                            codes_set(gid, 'time', t_m2dt.hour*100)
                            codes_set(gid, 'date', int(t_m2dt.strftime('%Y%m%d')))

                        codes_write(gid, f_handle)

                        # squeeze out information of last two steps
                        # contained in deac_vals[parId]
                        # Note that deac_vals[parId][0] has not been popped
                        # in this case

                        if step == c.maxstep and c.purefc or \
                           t_dt == t_enddate:
                            # last step
                            if c.purefc:
                                values = deac_vals[parId][3]
                                codes_set_values(gid, values)
                                codes_set(gid, 'stepRange', step)
                                #truedatetime = t_m2dt + timedelta(hours=2*int(c.dtime))
                                codes_write(gid, h_handle)
                            else:
                                codes_set(gid, 'stepRange', 0)
                                truedatetime = t_m2dt + timedelta(hours=2*int(c.dtime))
                                codes_set(gid, 'time', truedatetime.hour * 100)
                                codes_set(gid, 'date', int(truedatetime.strftime('%Y%m%d')))

                            if parId == 142 or parId == 143:
                                values = disaggregation.darain(list(reversed(deac_vals[parId])))
                            else:
                                values = disaggregation.dapoly(list(reversed(deac_vals[parId])))

                            # step before last step
                                codes_set(gid, 'stepRange', step-int(c.dtime))
                                #truedatetime = t_m2dt + timedelta(hours=int(c.dtime))
                                codes_write(gid, g_handle)
                            else:
                                truedatetime = t_m2dt + timedelta(hours=int(c.dtime))

                codes_release(gid)

                gid = codes_new_from_index(iid)

            f_handle.close()
            g_handle.close()
            h_handle.close()

        codes_index_release(iid)

            self._create_rr_grib_dummy(inputfiles.files[0], c.inputdir)

            self._prep_new_rrint(dims[0], dims[1], dims[2], lsp_np,
                                 cp_np, maxnum, index_keys, index_vals, c)

        return

    def _prep_new_rrint(self, ni, nj, nt, lsp_np, cp_np, maxnum, index_keys, index_vals, c):
        '''Calculates and writes out the disaggregated precipitation fields.

        Disaggregation is done in time and original times are written to the
        flux files, while the additional subgrid times are written to
        extra files output files. They are named like the original files with
        suffixes "_1" and "_2" for the first and second subgrid point.

        ni : int
            Amount of zonal grid points.

        nj : int
            Amount of meridional grid points.

        nt : int
            Number of time steps.

        lsp_np : numpy array of float
            The large scale precipitation fields for each time step.
            Shape (ni * nj, nt).

        cp_np : numpy array of float
            The convective precipitation fields for each time step.

        maxnum : int
            The maximum number of ensemble members. It is None
            if there are no or just one ensemble.

        index_keys : dictionary

        Return
        ------

        '''

        print('... disaggregation of precipitation with new method.')

        tmpfile = os.path.join(c.inputdir, 'rr_grib_dummy.grb')

        # initialize new numpy arrays for disaggregated fields
        if maxnum:
            lsp_new_np = np.zeros((maxnum, ni * nj, nt * 3), dtype=np.float64)
            cp_new_np = np.zeros((maxnum, ni * nj, nt * 3), dtype=np.float64)
        else:
            lsp_new_np = np.zeros((1, ni * nj, nt * 3), dtype=np.float64)
            cp_new_np = np.zeros((1, ni * nj, nt * 3), dtype=np.float64)

        # do the disaggregation, but neglect the last value of the
        # original time series. This one corresponds for example to
        # 24 hour, which we don't need. we use 0 - 23 UTC for a day.
        if maxnum:
            for inum in range(maxnum):
                for ix in range(ni*nj):
                    lsp_new_np[inum, ix, :] = disaggregation.IA3(lsp_np[inum, ix, :])[:-1]
                    cp_new_np[inum, ix, :] = disaggregation.IA3(cp_np[inum, ix, :])[:-1]
        else:
            for ix in range(ni*nj):
                lsp_new_np[0, ix, :] = disaggregation.IA3(lsp_np[ix, :])[:-1]
                cp_new_np[0, ix, :] = disaggregation.IA3(cp_np[ix, :])[:-1]

        # write to grib files (full/orig times to flux file and inbetween
        # times with step 1 and 2, respectively)
        print('... write disaggregated precipitation to files.')

        if maxnum:
            # empty set to save unique ensemble numbers which were already processed
            inumb = 0
        else:
            inumb = 0

        # index variable of disaggregated fields
        it = 0

            # or   prod = ('0'   , '20170505', '0', '12')
            #             (number, date      ,time, step)

            cdate = prod[index_keys.index('date')]
            ctime = '{:0>2}'.format(int(prod[index_keys.index('time')])//100)
            cstep = '{:0>3}'.format(int(prod[index_keys.index('step')]))

            date = datetime.strptime(cdate + ctime, '%Y%m%d%H')
            date += timedelta(hours=int(cstep))

            start_period, end_period = generate_retrieval_period_boundary(c)
            # skip all temporary times
            # which are outside the retrieval period
            if date < start_period or \
               date > end_period:
                continue

                inumb = int(prod[index_number])
                it = 0

            if maxnum:
            else:

            # per original time stamp: write original time step and
            # the two newly generated sub time steps
                fluxfilename = 'flux' + date.strftime('%Y%m%d.%H') + '.' + cstep
            else:
                fluxfilename = 'flux' + date.strftime('%Y%m%d%H') + numbersuffix

            # write original time step to flux file as usual
            fluxfile = GribUtil(os.path.join(c.inputdir, fluxfilename))
            fluxfile.set_keys(tmpfile, filemode='ab',
                              wherekeynames=['paramId'], wherekeyvalues=[142],
                              keynames=['perturbationNumber', 'date', 'time',
                                        'stepRange', 'values'],
                              keyvalues=[inumb, int(date.strftime('%Y%m%d')),
                                         date.hour*100, 0, lsp_new_np[inumb, :, it]],
                             )
                              wherekeynames=['paramId'], wherekeyvalues=[143],
                                         date.hour*100, 0, cp_new_np[inumb, :, it]]
                             )

            # rr for first subgrid point is identified by step = 1
                                         date.hour*100, '1', lsp_new_np[inumb, :, it+1]]
                             )
                                         date.hour*100, '1', cp_new_np[inumb, :, it+1]]
                             )

            # rr for second subgrid point is identified by step = 2
                                         date.hour*100, '2', lsp_new_np[inumb, :, it+2]]
                             )
                                         date.hour*100, '2', cp_new_np[inumb, :, it+2]]
                             )

            it = it + 3 # jump to next original time step in rr fields
        return

    def _create_rr_grib_dummy(self, ifile, inputdir):
        '''Creates a grib file with a dummy message for the two precipitation
        types lsp and cp each.

        ifile : str
            Filename of the input file to read the grib messages from.

            Path to the directory where the retrieved data is stored.

        Return
        ------

        '''

        gribfile = GribUtil(os.path.join(inputdir, 'rr_grib_dummy.grb'))

        gribfile.copy_dummy_msg(ifile, keynames=['paramId'],
                                keyvalues=[142], filemode='wb')

                                keyvalues=[143], filemode='ab')

        return

    def create(self, inputfiles, c):
        '''An index file will be created which depends on the combination
        of "date", "time" and "stepRange" values. This is used to iterate
        over all messages in each grib file which were passed through the
        parameter "inputfiles" to seperate specific parameters into fort.*
        files. Afterwards the FORTRAN program is called to convert
        the data fields all to the same grid and put them in one file
        per unique time step (combination of "date", "time" and
        "stepRange").

        Note
        ----
        This method is based on the ECMWF example index.py
        https://software.ecmwf.int/wiki/display/GRIB/index.py

        Return
        ------

        '''

        # generate start and end timestamp of the retrieval period
        start_period = datetime.strptime(c.start_date + c.time[0], '%Y%m%d%H')
        start_period = start_period + timedelta(hours=int(c.step[0]))
        end_period = datetime.strptime(c.end_date + c.time[-1], '%Y%m%d%H')
        end_period = end_period + timedelta(hours=int(c.step[-1]))

        # @WRF
        # THIS IS NOT YET CORRECTLY IMPLEMENTED !!!
        #
        # UNDER CONSTRUCTION !!!
        #
        #if c.wrf:
        #    table128 = init128(_config.PATH_GRIBTABLE)
        #    wrfpars = to_param_id('sp/mslp/skt/2t/10u/10v/2d/z/lsm/sst/ci/sd/\
        #                           stl1/stl2/stl3/stl4/swvl1/swvl2/swvl3/swvl4',
        #                          table128)

        # these numbers are indices for the temporary files "fort.xx"
        # which are used to seperate the grib fields to,
        # for the Fortran program input
        # 10: U,V | 11: T | 12: lnsp | 13: D | 16: sfc fields
        # 17: Q | 18: Q, SL, GG| 19: omega | 21: etadot | 22: clwc+ciwc
        fdict = {'10':None, '11':None, '12':None, '13':None, '16':None,
                 '17':None, '18':None, '19':None, '21':None, '22':None}

        iid = None

        # of grib messages via product
        else:
        # index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange

            #             (  date    ,time, step)

            print('current product: ', prod)

            if not gid:
                continue
#============================================================================================
            # remove old fort.* files and open new ones
            # they are just valid for a single product
            for k, f in fdict.items():
                fortfile = os.path.join(c.inputdir, 'fort.' + k)
                silent_remove(fortfile)
                fdict[k] = open(fortfile, 'wb')
            ctime = '{:0>2}'.format(codes_get(gid, 'time') // 100)
            cstep = '{:0>3}'.format(codes_get(gid, 'step'))
            timestamp = datetime.strptime(cdate + ctime, '%Y%m%d%H')
            timestamp += timedelta(hours=int(cstep))
            cdate_hour = datetime.strftime(timestamp, '%Y%m%d%H')

            # if basetime is used, adapt start/end date period
                time_delta = timedelta(hours=12-int(c.dtime))
                start_period = datetime.strptime(c.end_date + str(c.basetime),
                                               '%Y%m%d%H') - time_delta
                end_period = datetime.strptime(c.end_date + str(c.basetime),
                                             '%Y%m%d%H')

            if timestamp < start_period or \
               timestamp > end_period:
                continue

            # @WRF
            #
            #
            #if c.wrf:
            #    if 'olddate' not in locals() or cdate != olddate:
            #        fwrf = open(os.path.join(c.outputdir,
            #                    'WRF' + cdate + '.' + ctime + '.000.grb2'), 'wb')
            #        olddate = cdate[:]
            # savedfields remembers which fields were already used.
            savedfields = []
            # sum of cloud liquid and ice water content
            scwc = None
            while 1:
                if not gid:
                    break
                paramId = codes_get(gid, 'paramId')
                gridtype = codes_get(gid, 'gridType')
                if paramId == 77: # ETADOT
                    codes_write(gid, fdict['21'])
                elif paramId == 130: # T
                    codes_write(gid, fdict['11'])
                elif paramId == 131 or paramId == 132: # U, V wind component
                    codes_write(gid, fdict['10'])
                elif paramId == 133 and gridtype != 'reduced_gg': # Q
                    codes_write(gid, fdict['17'])
                elif paramId == 133 and gridtype == 'reduced_gg': # Q, gaussian
                    codes_write(gid, fdict['18'])
                elif paramId == 135: # W
                    codes_write(gid, fdict['19'])
                elif paramId == 152: # LNSP
                    codes_write(gid, fdict['12'])
                elif paramId == 155 and gridtype == 'sh': # D
                    codes_write(gid, fdict['13'])
                elif paramId == 246 or paramId == 247: # CLWC, CIWC
                    # sum cloud liquid water and ice
                    if scwc is None:
                        scwc = codes_get_values(gid)
                    else:
                        scwc += codes_get_values(gid)
                        codes_set_values(gid, scwc)
                        codes_set(gid, 'paramId', 201031)
                        codes_write(gid, fdict['22'])
                        scwc = None
                # @WRF
                # THIS IS NOT YET CORRECTLY IMPLEMENTED !!!
                #
                # UNDER CONSTRUCTION !!!
                #
                #elif c.wrf and paramId in [129, 138, 155] and \
                #      levtype == 'hybrid': # Z, VO, D
                #    # do not do anything right now
                #    # these are specific parameter for WRF
                #    pass
                else:
                    if paramId not in savedfields:
                        # SD/MSL/TCC/10U/10V/2T/2D/Z/LSM/SDOR/CVL/CVH/SR
                        # and all ADDPAR parameter
                        codes_write(gid, fdict['16'])
                        savedfields.append(paramId)
                    else:
                        print('duplicate ' + str(paramId) + ' not written')
                # @WRF
                #
                #
                #try:
                #    if c.wrf:
                #        # model layer
                #        if levtype == 'hybrid' and \
                #           paramId in [129, 130, 131, 132, 133, 138, 155]:
                #            codes_write(gid, fwrf)
                #        # sfc layer
                #        elif paramId in wrfpars:
                #except AttributeError:
                #    pass

            for f in fdict.values():
                f.close()
            # call for Fortran program to convert e.g. reduced_gg grids to
            # regular_ll and calculate detadot/dp
            pwd = os.getcwd()
            os.chdir(c.inputdir)
            if os.stat('fort.21').st_size == 0 and c.eta:
                print('Parameter 77 (etadot) is missing, most likely it is '
                      'not available for this type or date / time\n')
                print('Check parameters CLASS, TYPE, STREAM, START_DATE\n')
                my_error('fort.21 is empty while parameter eta '
                         'is set to 1 in CONTROL file')
# ============================================================================================
            # write out all output to log file before starting fortran programm
            sys.stdout.flush()

            # Fortran program creates file fort.15 (with u,v,etadot,t,sp,q)
            execute_subprocess([os.path.join(c.exedir,
                                             _config.FORTRAN_EXECUTABLE)],
                               error_msg='FORTRAN PROGRAM FAILED!')#shell=True)

            os.chdir(pwd)
            # create name of final output file, e.g. EN13040500 (ENYYMMDDHH)
                suffix = cdate[2:8] + '.' + ctime + '.' + cstep
            else:
                suffix = cdate_hour[2:10]

            if 'number' in index_keys:
                index_number = index_keys.index('number')
                if len(index_vals[index_number]) > 1:
                    suffix = suffix + '.N{:0>3}'.format(int(prod[index_number]))

            fnout = os.path.join(c.inputdir, c.prefix + suffix)
            # collect for final processing
            self.outputfilelist.append(os.path.basename(fnout))
            # # get additional precipitation subgrid data if available
            # if c.rrint:
                # self.outputfilelist.append(os.path.basename(fnout + '_1'))
                # self.outputfilelist.append(os.path.basename(fnout + '_2'))
            # create outputfile and copy all data from intermediate files
            # to the outputfile (final GRIB input files for FLEXPART)
            orolsm = os.path.basename(glob.glob(c.inputdir +
                                                '/OG_OROLSM__SL.*.' +
                                                c.ppid +
                                                '*')[0])
            fluxfile = 'flux' + cdate[0:2] + suffix
            if not c.cwc:
                flist = ['fort.15', fluxfile, 'fort.16', orolsm]
            else:
                flist = ['fort.15', 'fort.22', fluxfile, 'fort.16', orolsm]

            with open(fnout, 'wb') as fout:
                for f in flist:
                    shutil.copyfileobj(open(os.path.join(c.inputdir, f), 'rb'),
                                       fout)

            if c.omega:
                with open(os.path.join(c.outputdir, 'OMEGA'), 'wb') as fout:
                    shutil.copyfileobj(open(os.path.join(c.inputdir, 'fort.25'),
                                            'rb'), fout)

        # @WRF
        #
        #
        #if c.wrf:
        #    fwrf.close()

        return

    def calc_extra_elda(self, path, prefix):
        ''' Calculates extra ensemble members for ELDA - Stream.

        This is a specific feature which doubles the number of ensemble members
        for the ELDA Stream.

        path : str
            Path to the output files.

        prefix : str
            The prefix of the output filenames as defined in Control file.

        Return
        ------

        '''
        from eccodes import (codes_grib_new_from_file, codes_get_array,
                             codes_set_array, codes_release,
                             codes_set, codes_write)

        # max number
        maxnum = int(self.number.split('/')[-1])

        # get a list of all prepared output files with control forecast (CF)
        cf_filelist = UioFiles(path, prefix + '*.N000')
        cf_filelist.files = sorted(cf_filelist.files)

        for cffile in cf_filelist.files:
            with open(cffile, 'rb') as f:
                cfvalues = []
                while True:
                    fid = codes_grib_new_from_file(f)
                    if fid is None:
                        break
                    cfvalues.append(codes_get_array(fid, 'values'))
                    codes_release(fid)

            filename = cffile.split('N000')[0]
            for i in range(1, maxnum + 1):
                # read an ensemble member
                g = open(filename + 'N{:0>3}'.format(i), 'rb')
                # create file for newly calculated ensemble member
                h = open(filename + 'N{:0>3}'.format(i+maxnum), 'wb')
                # number of message in grib file
                j = 0
                    gid = codes_grib_new_from_file(g)
                    if gid is None:
                        break
                    values = codes_get_array(gid, 'values')
                    # generate a new ensemble member by subtracting
                    # 2 * ( current time step value - last time step value )
                    codes_set_array(gid, 'values',
                                    values-2*(values-cfvalues[j]))
                    codes_set(gid, 'number', i+maxnum)
                    codes_write(gid, h)
                    codes_release(gid)
                    j += 1

                g.close()
                h.close()
                print('wrote ' + filename + 'N{:0>3}'.format(i+maxnum))
                self.outputfilelist.append(
                    os.path.basename(filename + 'N{:0>3}'.format(i+maxnum)))

        return

    def process_output(self, c):
        '''Postprocessing of FLEXPART input files.

        The grib files are postprocessed depending on the selection in
        CONTROL file. The resulting files are moved to the output
        directory if its not equal to the input directory.
        The following modifications might be done if
        properly switched in CONTROL file:
        GRIB2 - Conversion to GRIB2
        ECTRANS - Transfer of files to gateway server
        ECSTORAGE - Storage at ECMWF server

        Return
        ------

        '''

        print('\n\nPostprocessing:\n Format: {}\n'.format(c.format))

        if _config.FLAG_ON_ECMWFSERVER:
            print('ecstorage: {}\n ecfsdir: {}\n'.
                  format(c.ecstorage, c.ecfsdir))
            print('ectrans: {}\n gateway: {}\n destination: {}\n '
                  .format(c.ectrans, c.gateway, c.destination))

        print('Output filelist: ')
        print(sorted(self.outputfilelist))

        for ofile in self.outputfilelist:
            ofile = os.path.join(self.inputdir, ofile)

            if c.format.lower() == 'grib2':
                execute_subprocess(['grib_set', '-s', 'edition=2,' +
                                    'productDefinitionTemplateNumber=8',
                                    ofile, ofile + '_2'],
                                   error_msg='GRIB2 CONVERSION FAILED!')

                execute_subprocess(['mv', ofile + '_2', ofile],
                                   error_msg='RENAMING FOR NEW GRIB2 FORMAT '
                                   'FILES FAILED!')

            if c.ectrans and _config.FLAG_ON_ECMWFSERVER:
                execute_subprocess(['ectrans', '-overwrite', '-gateway',
                                    c.gateway, '-remote', c.destination,
                                    '-source', ofile],
                                   error_msg='TRANSFER TO LOCAL SERVER FAILED!')

            if c.ecstorage and _config.FLAG_ON_ECMWFSERVER:
                execute_subprocess(['ecp', '-o', ofile,
                                    os.path.expandvars(c.ecfsdir)],
                                   error_msg='COPY OF FILES TO ECSTORAGE '
                                   'AREA FAILED!')

            if c.outputdir != c.inputdir:
                execute_subprocess(['mv', os.path.join(c.inputdir, ofile),
                                    c.outputdir],
                                   error_msg='RELOCATION OF OUTPUT FILES '
                                   'TO OUTPUTDIR FAILED!')

        return

#
# @Date: July 2014
#
#        - changed some naming
#
#
#

# MODULES

import os

# CLASS
class GribUtil(object):
    '''
    Class for GRIB utilities (new methods) based on GRIB API

    The GRIB API provides all necessary tools to work directly with the
    grib files. Nevertheless, the GRIB API tools are very basic and are in
    direct connection with the grib files. This class provides some higher
    functions which apply a set of GRIB API tools together in the respective
    context. So, the class initially contains a list of grib files (their
    names) and the using program then applies the methods directly on the
    class objects without having to think about how the actual GRIB API
    tools have to be arranged.
    '''
    def __init__(self, filenames):
        '''Initialise an object of GribUtil and assign a list of filenames.

        filenames : :obj:`list` of :obj:`strings`
             A list of filenames.

        Return
        ------

        '''

        self.filenames = filenames

        return

    def get_keys(self, keynames, wherekeynames, wherekeyvalues):
        '''Get keyvalues for a given list of keynames a where statement
        can be given (list of key and list of values)

        keynames : :obj:`list` of :obj:`string`
            List of keynames.

        wherekeynames : :obj:`list` of :obj:`string`
            List of key names for indexing grib message parameter.

        wherekeyvalues : :obj:`list` of :obj:`string`
            List of key values corresponding the key names.

        Return
        ------
        return_list : :obj:`list` of :obj:`string`
            List of keyvalues for given keynames.
        '''
        from eccodes import (codes_new_from_file, codes_is_defined, codes_get,
                             codes_release)

        return_list = []

        with open(self.filenames, 'rb') as fileid:

            gid = codes_new_from_file(fileid)

            if len(wherekeynames) != len(wherekeyvalues):
                raise Exception("Number of key values and key names must be \
                                 the same. Give a value for each keyname!")

            select = True
            i = 0
            for wherekey in wherekeynames:
                if not codes_is_defined(gid, wherekey):
                    raise Exception("where key was not defined")

                select = (select and (str(wherekeyvalues[i]) ==
                                      str(codes_get(gid, wherekey))))
                i += 1

            if select:
                llist = []
                for key in keynames:
                    llist.extend([str(codes_get(gid, key))])
                return_list.append(llist)

            codes_release(gid)

        return return_list

    def set_keys(self, fromfile, keynames, keyvalues, wherekeynames,
                 wherekeyvalues, filemode='wb'):
        '''Opens the file to read the grib messages and then write
        the selected messages (with wherekeys) to a new output file.
        Also, the keyvalues of the passed list of keynames are set.

        fromfile : :obj:`string`

            List of keynames to set in the selected messages.
            Default is an empty list.

        keyvalues : :obj:`list` of :obj:`string`
            List of keyvalues to set in the selected messages.

            List of keynames to select correct message.

            List of keyvalues for keynames to select correct message.

        filemode : :obj:`string`, optional
            Sets the mode for the output file. Default is "wb".

        Return
        ------

        '''
        from eccodes import (codes_grib_new_from_file, codes_is_defined,
                             codes_get, codes_set, codes_write,
                             codes_set_values, codes_release)

        if len(wherekeynames) != len(wherekeyvalues):
            raise Exception("Give a value for each keyname!")

        fout = open(self.filenames, filemode)

        with open(fromfile, 'rb') as fin:
            gid = codes_grib_new_from_file(fin)

            i = 0
                    raise Exception("wherekey was not defined")

                i += 1

            if select:
                i = 0
                    if key == 'values':
                        codes_set_values(gid, keyvalues[i])
                    else:
                        codes_set(gid, key, keyvalues[i])
                    i += 1

                codes_write(gid, fout)

        fout.close()

        return

    def copy_dummy_msg(self, filename_in, keynames, keyvalues,
                       selectwhere=True, filemode='wb'):
        '''Add the content of another input grib file to the objects file but
        only messages corresponding to keys/values passed to the function.
        The selectWhere switch decides if to copy the keys equal to (True) or
        different to (False) the keynames/keyvalues list passed to the function.

        filename_in : :obj:`string`

        selectwhere : :obj:`boolean`, optional
            Decides if to copy the keynames and values equal to (True) or
            different to (False) the keynames/keyvalues list passed to the
            function. Default is True.

            List of keyvalues.

        Return
        ------

        '''
                             codes_get, codes_release, codes_write)

        if len(keynames) != len(keyvalues):

        fields = 0

        with open(filename_in, 'rb') as fin:
            if fields >= 1:
                fout.close()
                return

            i = 0
            for key in keynames:
                if not codes_is_defined(gid, key):
                    raise Exception("Key was not defined")

                if selectwhere:
                    select = (select and (str(keyvalues[i]) ==
                                          str(codes_get(gid, key))))
                else:
                    select = (select and (str(keyvalues[i]) !=
                i += 1

            if select:
                fields = fields + 1

        return

    def index(self, index_keys, index_file="my.idx"):
        '''Create index file from a list of files if it does not exist or
        read an index file.

        index_keys: :obj:`list` of :obj:`string`
            Contains the list of key parameter names from
            which the index is to be created.

        index_file: :obj:`string`, optional
            Filename where the indices are stored.
            Default is "my.idx".

        Return
        ------
        iid: :obj:`integer`
            Grib index id.
        '''
        from eccodes import (codes_index_read, codes_index_new_from_file,
                             codes_index_add_file, codes_index_write)

        print("... index will be done")
        iid = None

        if os.path.exists(index_file):
            iid = codes_index_read(index_file)
            print("Use existing index file: %s " % (index_file))
        else:
            for filename in self.filenames:
                print("Inputfile: %s " % (filename))
                if iid is None:
                    iid = codes_index_new_from_file(filename, index_keys)
                else:
                    codes_index_add_file(iid, filename)

            if iid is not None:
                codes_index_write(iid, index_file)

        print('... index done')

        return iid

#
#
#
#   November 2015 - Leopold Haimberger (University of Vienna):
#        - optimized display_info
#        - optimized data_retrieve and seperate between python and shell
#          script call
#
#        - added writing of mars request attributes to a csv file
#
#
#

# MODULES

import os
import sys
import traceback

#pylint: disable=invalid-name
try:
    ec_api = True
    import ecmwfapi
except ImportError:
    ec_api = False

try:
    cds_api = True
    import cdsapi
    cds_api = False
#pylint: enable=invalid-name
# CLASS
class MarsRetrieval(object):
    '''Specific syntax and content for submission of MARS retrievals.

    A MARS revtrieval has a specific syntax with a selection of keywords and
    their corresponding values. This class provides the necessary functions
    by displaying the selected parameters and their values and the actual
    retrievement of the data through a mars request or a Python web api
    interface. The initialization already expects all the keyword values.

    A description of MARS keywords/arguments and examples of their
    values can be found here:
    https://software.ecmwf.int/wiki/display/UDOC/\
                   Identification+keywords#Identificationkeywords-class

    marsclass : str, optional

    dataset : str, optional
        data.

    type : str, optional
        Determines the type of fields to be retrieved.

    levtype : str, optional
        Denotes type of level.

    levelist : str, optional

    repres : str, optional
        Selects the representation of the archived data.

    date : str, optional
        Specifies the Analysis date, the Forecast base date or
        Observations date.

    resol : str, optional

    stream : str, optional

    area : str, optional

    time : str, optional
        Specifies the time of the data in hours and minutes.

    step : str, optional
        Specifies the forecast time step from forecast base time.

    expver : str, optional

    number : str, optional

    accuracy : str, optional

    grid : str, optional

    gaussian : str, optional

    target : str, optional
        Specifies a file into which data is to be written after
        retrieval or manipulation.

    param : str, optional
        Specifies the meteorological parameter.
    '''

    def __init__(self, server, public, marsclass="EA", dataset="", type="",
                 levtype="", levelist="", repres="", date="", resol="",
                 stream="", area="", time="", step="", expver="1",
                 number="", accuracy="", grid="", gaussian="", target="",
                 param=""):
        '''Initialises the instance of the MarsRetrieval class and
        defines and assigns a set of the necessary retrieval parameters
        for the FLEXPART input data.
        A description of MARS keywords/arguments, their dependencies
        on each other and examples of their values can be found here:

        https://software.ecmwf.int/wiki/display/UDOC/MARS+keywords

            This is the connection to the ECMWF data servers.
            It is needed for the pythonic access of ECMWF data.

        public : int
            Decides which Web API version is used:
            0: member-state users and full archive access
            1: public access and limited access to the public server and
               datasets. Needs the parameter dataset.
            Default is "0" and for member-state users.

        marsclass : str, optional
            Characterisation of dataset. E.g. EI (ERA-Interim),
            E4 (ERA40), OD (Operational archive), EA (ERA5).
            Default is the ERA5 dataset "EA".

        dataset : str, optional
            For public datasets there is the specific naming and parameter
            dataset which has to be used to characterize the type of
            data. Usually there is less data available, either in times,
            domain or parameter.
            Default is an empty string.

        type : str, optional
            Determines the type of fields to be retrieved.
            Selects between observations, images or fields.
            Examples for fields: Analysis (an), Forecast (fc),
            Perturbed Forecast (pf), Control Forecast (cf) and so on.

        levtype : str, optional
            Denotes type of level. Has a direct implication on valid
            levelist values!
            E.g. model level (ml), pressure level (pl), surface (sfc),
            potential vorticity (pv), potential temperature (pt)
            and depth (dp).

        levelist : str, optional
            Specifies the required levels. It has to have a valid
            correspondence to the selected levtype.
            Examples: model level: 1/to/137, pressure levels: 500/to/1000

        repres : str, optional
            Selects the representation of the archived data.
            E.g. sh - spherical harmonics, gg - Gaussian grid,
            ll - latitude/longitude, ...

        date : str, optional
            Specifies the Analysis date, the Forecast base date or
            Observations date. Valid formats are:
            Absolute as YYYY-MM-DD or YYYYMMDD.

        resol : str, optional
            Specifies the desired triangular truncation of retrieved data,
            before carrying out any other selected post-processing.
            The default is automatic truncation (auto), by which the lowest
            resolution compatible with the value specified in grid is
            automatically selected for the retrieval.
            Users wanting to perform post-processing from full spectral
            resolution should specify Archived Value (av).
            The following are examples of existing resolutions found in
            the archive: 63, 106, 159, 213, 255, 319, 399, 511, 799 or 1279.
            This keyword has no meaning/effect if the archived data is
            not in spherical harmonics representation.
            The best selection can be found here:
            https://software.ecmwf.int/wiki/display/UDOC/\
                  Retrieve#Retrieve-Truncationbeforeinterpolation

        stream : str, optional
            Identifies the forecasting system used to generate the data.
            E.g. oper (Atmospheric model), enfo (Ensemble forecats), ...

        area : str, optional
            Specifies the desired sub-area of data to be extracted.
            Areas can be defined to wrap around the globe.

            Latitude values must be given as signed numbers, with:
                north latitudes (i.e. north of the equator)
                    being positive (e.g: 40.5)
                south latitutes (i.e. south of the equator)
                    being negative (e.g: -50.5)
            Longtitude values must be given as signed numbers, with:
                east longitudes (i.e. east of the 0 degree meridian)
                    being positive (e.g: 35.0)
                west longitudes (i.e. west of the 0 degree meridian)
                    being negative (e.g: -20.5)

            E.g.: North/West/South/East

        time : str, optional
            Specifies the time of the data in hours and minutes.
            Valid values depend on the type of data: Analysis time,
            Forecast base time or First guess verification time
            (all usually at synoptic hours: 00, 06, 12 and 18 ).
            Observation time (any combination in hours and minutes is valid,
            subject to data availability in the archive).
            The syntax is HHMM or HH:MM. If MM is omitted it defaults to 00.

        step : str, optional
            Valid values are hours (HH) from forecast base time. It also
            specifies the length of the forecast which verifies at
            First Guess time.
            E.g. 1/3/6-hourly

        expver : str, optional
            The version of the dataset. Each experiment is assigned a
            unique code (version). Production data is assigned 1 or 2,
            and experimental data in Operations 11, 12 ,...
            Research or Member State's experiments have a four letter
            experiment identifier.
            Default is "1".

        number : str, optional
            Selects the member in ensemble forecast run. (Only then it
            is necessary.) It has a different meaning depending on
            the type of data.
            E.g. Perturbed Forecasts: specifies the Ensemble forecast member

        accuracy : str, optional
            Specifies the number of bits per value to be used in the
            generated GRIB coded fields.
            A positive integer may be given to specify the preferred number
            of bits per packed value. This must not be greater than the
            number of bits normally used for a Fortran integer on the
            processor handling the request (typically 32 or 64 bit).
            Within a compute request the accuracy of the original fields
            can be passed to the result field by specifying accuracy=av.

        grid : str, optional
            Specifies the output grid which can be either a Gaussian grid
            or a Latitude/Longitude grid. MARS requests specifying
            grid=av will return the archived model grid.

            Lat/Lon grid: The grid spacing needs to be an integer
            fraction of 90 degrees e.g. grid = 0.5/0.5

            Gaussian grid: specified by a letter denoting the type of
            Gaussian grid followed by an integer (the grid number)
            representing the number of lines between the Pole and Equator,
            e.g.
            grid = F160 - full (or regular) Gaussian grid with
                   160 latitude lines between the pole and equator
            grid = N320 - ECMWF original reduced Gaussian grid with
                   320 latitude lines between the pole and equator,
                   see Reduced Gaussian Grids for grid numbers used at ECMWF
            grid = O640 - ECMWF octahedral (reduced) Gaussian grid with
                   640 latitude lines between the pole and equator

        gaussian : str, optional
            This parameter is deprecated and should no longer be used.
            Specifies the desired type of Gaussian grid for the output.
            Valid Gaussian grids are quasi-regular (reduced) or regular.
            Keyword gaussian can only be specified together with
            keyword grid. Gaussian without grid has no effect.

        target : str, optional
            Specifies a file into which data is to be written after
            retrieval or manipulation. Path names should always be
            enclosed in double quotes. The MARS client supports automatic
            generation of multiple target files using MARS keywords
            enclosed in square brackets [ ].  If the environment variable
            MARS_MULTITARGET_STRICT_FORMAT is set to 1 before calling mars,
            the keyword values will be used in the filename as shown by
            the ecCodes GRIB tool grib_ls -m, e.g. with
            MARS_MULTITARGET_STRICT_FORMAT set to 1 the keywords time,
            expver and param will be formatted as 0600, 0001 and 129.128
            rather than 600, 1 and 129.

        param : str, optional
            Specifies the meteorological parameter.
            The list of meteorological parameters in MARS is extensive.
            Their availability is directly related to their meteorological
            meaning and, therefore, the rest of directives specified
            in the MARS request.
            Meteorological parameters can be specified by their
            GRIB code (param=130), their mnemonic (param=t) or
            full name (param=temperature).
            The list of parameter should be seperated by a "/"-sign.
            E.g. 130/131/133

        Return
        ------

        '''

        self.marsclass = marsclass
        self.dataset = dataset
        self.type = type
        self.levtype = levtype
        self.levelist = levelist
        self.repres = repres
        self.date = date
        self.resol = resol
        self.stream = stream
        self.area = area
        self.time = time
        self.step = step
        self.expver = expver
        self.number = number
        self.accuracy = accuracy
        self.grid = grid
        self.gaussian = gaussian
        self.target = target
        self.param = param

        return

    def display_info(self):
        '''Prints all class attributes and their values to the
        standard output.

        Return
        ------

        '''
        # Get all class attributes and their values as a dictionary

        # iterate through all attributes and print them
        # with their corresponding values
            if item[0] in ['server', 'public']:
                pass
            else:
                print(item[0] + ': ' + str(item[1]))

        return

    def print_infodata_csv(self, inputdir, request_number):
        '''Write all request parameter in alpabetical order into a "csv" file.

        inputdir : str
            The path where all data from the retrievals are stored.

        request_number : int
            Number of mars requests for flux and non-flux data.

        Return
        ------

        '''

        del attrs['server']
        del attrs['public']

        # open a file to store all requests to
        with open(os.path.join(inputdir,
                               _config.FILE_MARS_REQUESTS), 'a') as f:
            f.write(str(request_number) + ', ')
            f.write(', '.join(str(attrs[key])
                              for key in sorted(attrs.keys())))
            f.write('\n')

        return

    def data_retrieve(self):
        '''Submits a MARS retrieval. Depending on the existence of
        ECMWF Web-API or CDS API it is submitted via Python or a
        subprocess in the Shell. The parameter for the mars retrieval
        are taken from the defined class attributes.

        Return
        ------

        '''

        # eliminate unnecessary attributes from the dictionary attrs

        # exchange parameter name for marsclass
        mclass = attrs.get('marsclass')
        del attrs['marsclass']
        attrs['class'] = mclass

        # prepare target variable as needed for the Web API or CDS API mode
        # within the dictionary for full access
        # as a single variable for public access
        target = attrs.get('target')
        if not int(self.public):
            del attrs['target']
        print('target: ' + target)

        # find all keys without a value and convert all other values to strings
        empty_keys = []
        for key, value in attrs.items():
            if value == '':
                empty_keys.append(str(key))
            else:
                attrs[key] = str(value)

        # delete all empty parameter from the dictionary
        for key in empty_keys:
            del attrs[key]

#        attrs['ppengine'] = 'emos'

        # MARS request via Python script
        if self.server:
            try:
                if cds_api and isinstance(self.server, cdsapi.Client):
                    print('RETRIEVE ERA5 WITH CDS API!')
                    self.server.retrieve(_config.CDS_DATASET,
                                         attrs, target)
                elif ec_api and isinstance(self.server, ecmwfapi.ECMWFDataServer):
                    print('RETRIEVE PUBLIC DATA (NOT ERA5)!')
                    self.server.retrieve(attrs)
                elif ec_api and isinstance(self.server, ecmwfapi.ECMWFService):
                    print('EXECUTE NON-PUBLIC RETRIEVAL (NOT ERA5)!')
                    self.server.execute(attrs, target)
                else:
                    print('ERROR:')
                    print('No match for Web API instance!')
                    raise IOError
            except Exception as e:
                print('\n\nMARS Request failed!')
                print(e)
                print(traceback.format_exc())
                sys.exit()

        # MARS request via call in shell
        else:
            request_str = 'ret'
            for key, value in attrs.items():
                request_str = request_str + ',' + key + '=' + str(value)
            request_str += ',target="' + target + '"'
            p = subprocess.Popen(['mars'], #'-e'],
                                 stdin=subprocess.PIPE,
                                 stdout=subprocess.PIPE,
                                 stderr=subprocess.PIPE,
                                 bufsize=1)
            pout = p.communicate(input=request_str.encode())[0]
            print(pout.decode())

            if 'Some errors reported' in pout.decode():
                print('MARS Request failed - please check request')
                raise IOError
            elif os.stat(target).st_size == 0:
                print('MARS Request returned no data - please check request')

        return

#
#
#
#        - modified method list_files to work with glob instead of listdir
#        - added pattern search in method list_files
#
#    February - December 2018 - Anne Philipp (University of Vienna):
#        - optimisation of method list_files since it didn't work correctly
#          for sub directories
#        - additional speed up of method list_files
#        - modified the class so that it is initiated with a pattern instead
#          of suffixes. Gives more precision in selection of files.
#        - added delete method
#
#
#

# MODULES
import os
import sys
import fnmatch

# software specific modules from flex_extract
from Mods.tools import silent_remove, get_list_as_string

# CLASS

class UioFiles(object):
    """Collection of files matching a specific pattern.

    The pattern can contain regular expressions for the files.
    The files are listed and can be transformed to a single string or
    they can be deleted.

    path : str
        Directory where to list the files.

    pattern : str
        Regular expression pattern. For example: '*.grb'

    files : list of str
        List of files matching the pattern in the path.
    """
    # CLASS METHODS
    def __init__(self, path, pattern):
        """Assignes a specific pattern for these files.

        path : str
            Directory where to list the files.

        pattern : str
            Regular expression pattern. For example: '*.grb'

        Return
        ------

        """

        self.path = path
        self.pattern = pattern
        self.files = []

        self._list_files(self.path)

        return

    def _list_files(self, path):
        """Lists all files in the directory with the matching
        regular expression pattern.

        path : str
            Path to the files.

        Return
        ------

        """
        # Get the absolute path
        path = os.path.abspath(path)

        # get all files in the dir and subdir as absolut path
        # pylint: disable=W0612
        for root, dirnames, filenames in os.walk(path):
            for filename in fnmatch.filter(filenames, self.pattern):
                self.files.append(os.path.join(root, filename))

        return

        """Converts the list of files into a single string.
        The entries are sepereated by "," sign.

        Return
        ------
        files_string : str
            The content of the list as a single string.
        """

        filenames = [os.path.basename(f) for f in self.files]
        files_string = get_list_as_string(filenames, concatenate_sign=', ')

        return files_string

    def delete_files(self):
        """Deletes the files.

        Return
        ------

        """

        for old_file in self.files:
            silent_remove(old_file)

        return

#
#
#
#        - moved the getEIdata program into a function "get_mars_data"
#        - moved the AgurmentParser into a separate function
#        - adapted the function for use in flex_extract
#        - renamed source file to get_mars_data
#
#        - added structured documentation
#        - minor changes in programming style for consistence
#        - added function main and moved function calls vom __main__ there
#          (necessary for better documentation with docstrings for later
#          online documentation)
#        - use of UIFiles class for file selection and deletion
#        - separated get_mars_data function into several smaller pieces:
#          write_reqheader, mk_server, mk_dates, remove_old, do_retrievment
#
#
#
'''This script extracts MARS data from ECMWF.

At first, the necessary parameters from command line and CONTROL files are
extracted. They define the data set to be extracted from MARS.

This file can also be imported as a module and contains the following

    * main            - the main function of the script
    * get_mars_data   - overall control of ECMWF data retrievment
    * write_reqheader - writes the header into the mars_request file
    * mk_server       - creates the server connection to ECMWF servers
    * mk_dates        - defines the start and end date
    * remove_old      - deletes old retrieved grib files
    * do_retrieval    - creates individual retrievals

Type get_mars_data.py --help
'''
# MODULES

import os
import sys
import inspect

# software-specific classes and modules from flex_extract
# add path to local main Python path for flex_extract to get full access
sys.path.append(os.path.dirname(os.path.abspath(
    inspect.getfile(inspect.currentframe()))) + '/../')
# pylint: disable=wrong-import-position
from Mods.tools import (setup_controldata, my_error, normal_exit, make_dir)
from Classes.EcFlexpart import EcFlexpart
# pylint: enable=wrong-import-position
# pylint: disable=invalid-name
try:

try:
# pylint: enable=invalid-name
# FUNCTION
    '''Controls the program to retrieve data from MARS.

    This is done if called directly from command line.
    Then, arguments and control file are taken as input.

    Return
    ------

    '''

    c, _, _, _ = setup_controldata()
    get_mars_data(c)
    normal_exit('Retrieving MARS data: Done!')

    return

def get_mars_data(c):
    '''Retrieves the ECMWF data required for a FLEXPART simulation.

    Start and end dates for retrieval period are set. Retrievals
    are divided into shorter periods if necessary and if datechunk parameter
    is set.

    Return
    ------

    '''
    c.ec_api = ec_api
    c.cds_api = cds_api

    if not os.path.exists(c.inputdir):
        make_dir(c.inputdir)

    if c.request == 0:
        print("Retrieving ECMWF data!")
    else:
        if c.request == 1:
            print("Printing MARS requests!")
        elif c.request == 2:
            print("Retrieving ECMWF data and printing MARS request!")
        write_reqheader(os.path.join(c.inputdir, _config.FILE_MARS_REQUESTS))

    print("start date %s " % (c.start_date))
    print("end date %s " % (c.end_date))

    server = mk_server(c)

    # if data are to be retrieved, clean up any old grib files
    if c.request == 0 or c.request == 2:
        remove_old('*grb', c.inputdir)

    # --------------  flux data ------------------------------------------------
    start, end, datechunk = mk_dates(c, fluxes=True)
    do_retrievement(c, server, start, end, datechunk, fluxes=True)

    # --------------  non flux data --------------------------------------------
    start, end, datechunk = mk_dates(c, fluxes=False)
    do_retrievement(c, server, start, end, datechunk, fluxes=False)

    return

def write_reqheader(marsfile):
    '''Writes header with column names into MARS request file.

    marsfile : str
        Path to the MARS request file.

    Return
    ------

    '''
    MR = MarsRetrieval(None, None)
    attrs = vars(MR).copy()
    del attrs['server']
    del attrs['public']
    with open(marsfile, 'w') as f:
        f.write('request_number' + ', ')
        f.write(', '.join(str(key) for key in sorted(attrs.keys())))
        f.write('\n')

    return

def mk_server(c):
    '''Creates a server connection with available Python API.

    The API selected depends on availability and the data set to be retrieved.
    The CDS API is used for ERA5 data, no matter whether the user is a 
    member-state or a public user. 
    ECMWF WebAPI is used for all other available datasets.

    Return
    ------
    server : ECMWFDataServer, ECMWFService or Client
        Connection to ECMWF server via python interface ECMWF WebAPI or CDS API.

    '''
    if cds_api and (c.marsclass.upper() == 'EA'):
        server = cdsapi.Client()
        c.ec_api = False
    elif c.ec_api:
        if c.public:
            server = ecmwfapi.ECMWFDataServer()
        else:
            server = ecmwfapi.ECMWFService("mars")
        c.cds_api = False
    else:
        server = False

    print('Using ECMWF WebAPI: ' + str(c.ec_api))
    print('Using CDS API: ' + str(c.cds_api))

    return server

def check_dates_for_nonflux_fc_times(types, times):
    '''Checks if the time 18UTC corresponds to forecast field.

    types : list of str
        List of field types.

    times : list of str or str
        The time in hours of the field.

    Return
    ------
    True or False

    '''
    for ty, ti in zip(types, times):
        if ty.upper() == 'FC' and int(ti) == 18:
            return True
    return False

def mk_dates(c, fluxes):
    '''Prepares start and end date depending on flux or non-flux type of data.

    If forecasts for a maximum of one day (24 h) are to be retrieved, then
    collect accumulation data (flux data) with additional days in the
    beginning and at the end (needed for complete disaggregation of
    original period)

    If forecast data for more than +24 h are to be retrieved, then
    collect accumulation data (flux data) with the exact start and end date
    (disaggregation will be done for the exact time period with
    boundary conditions)

    Since for basetime the extraction contains the 12 hours upfront,
    if basetime is 0, the starting date has to be the day before

    fluxes : boolean, optional
        Decides if the flux parameter settings are stored or
        the rest of the parameter list.
        Default value is False.

    Return
    ------
    start : datetime
        The start date of the retrieving data set.

    end : datetime
        The end date of the retrieving data set.

    chunk : datetime
        Time period in days for one single mars retrieval.

    '''
    start = datetime.strptime(c.start_date, '%Y%m%d')
    end = datetime.strptime(c.end_date, '%Y%m%d')
    chunk = timedelta(days=int(c.date_chunk))

    if c.basetime == 0:
        start = start - timedelta(days=1)

    if c.purefc and fluxes and c.maxstep < 24:
        end = end + timedelta(days=1)

    if not c.purefc and fluxes and not c.basetime == 0:

    # if we have non-flux forecast data starting at 18 UTC
    # we need to start retrieving data one day in advance
    if not fluxes and check_dates_for_nonflux_fc_times(c.type, c.time):

    return start, end, chunk

def remove_old(pattern, inputdir):
    '''Deletes old retrieval files from current input directory
    matching the pattern.

        The substring pattern which identifies the files to be deleted.

    inputdir : str, optional
        Path to the directory where the retrieved data are stored.

    Return
    ------

    '''
    print('... removing old files in ' + inputdir)

    tobecleaned = UioFiles(inputdir, pattern)
    tobecleaned.delete_files()

    return

def do_retrievement(c, server, start, end, delta_t, fluxes=False):
    '''Divides the total retrieval period into smaller chunks and
    retrieves the data from MARS.

            The server connection to ECMWF.

        The start date of the retrieval.

        The end date of the retrieval.

    delta_t : datetime
        Delta_t + 1 is the maximum time period of a single retrieval.

        Decides if the flux parameters are to be retrieved or

    Return
    ------

    '''

    # since actual day also counts as one day,
    # we only need to add datechunk - 1 days to retrieval for a period
    delta_t_m1 = delta_t - timedelta(days=1)

    day = start
    while day <= end:
        flexpart = EcFlexpart(c, fluxes)
        tmpday = day + delta_t_m1
        if tmpday < end:
            dates = day.strftime("%Y%m%d") + "/to/" + \
                    tmpday.strftime("%Y%m%d")
        else:
                    end.strftime("%Y%m%d")

        print("... retrieve " + dates + " in dir " + c.inputdir)

        try:
            flexpart.retrieve(server, dates, c.public, c.request, c.inputdir)
            my_error('MARS request failed')

        day += delta_t

    return

    main()

#
#
#
#        - using the WebAPI also for general MARS retrievals
#        - dividing retrievals of longer periods into digestable chunks
#        - retrieve also longer-term forecasts, not only analyses and
#          short-term forecast data
#        - conversion into GRIB2
#       
#
#        - BUGFIX: removed call of clean_up-Function after call of
#               prepareFlexpart in main since it is already called in
#               prepareFlexpart at the end!
#        - created function main and moved the two function calls for
#          arguments and prepare_flexpart into it
#
#
#
# *******************************************************************************
# pylint: disable=ungrouped-imports
'''This script prepares the final version of the grib files which are
then used by FLEXPART.

It converts the bunch of grib files extracted via get_mars_data before,
by doing the necessary conversion to obtain consistent grids or the
disaggregation of flux data. Finally, the data fields are combined
in files per hour available with the naming convention xxYYMMDDHH,
where xx should be 2 arbitrary letters (mostly xx is chosen to be "EN").

    * main
    * prepare_flexpart

Type: prepare_flexpart.py --help
'''

# MODULES

import datetime
import os
import sys

# add path to local main python path for flex_extract to get full access
#import _config
from Mods.checks import check_ppid
from Mods.tools import (setup_controldata, clean_up, make_dir, normal_exit)

# FUNCTION
    '''Controls the program to prepare FLEXPART input files from MARS data.

    This is done if called directly from the command line.

    Return
    ------

    '''

    c, ppid, _, _ = setup_controldata()
    prepare_flexpart(ppid, c)
    normal_exit('Preparing FLEXPART output files: Done!')

    return

def prepare_flexpart(ppid, c):
    '''Converts the MARS data into files ready as input for FLEXPART.

    Certain fields are converted to a different grid and the flux
    data are disaggregated. Fields are collected by hour and stored in a file 
    with a specific naming convention.

    ppid : int
        Contains the ppid number of the current ECMWF job. It will be None if
        the method was called within this module.

    Return
    ------

    '''
    check_ppid(c, ppid)

    # create the start and end date
    start = datetime.date(year=int(c.start_date[:4]),
                          month=int(c.start_date[4:6]),
                          day=int(c.start_date[6:]))

    end = datetime.date(year=int(c.end_date[:4]),
                        month=int(c.end_date[4:6]),
                        day=int(c.end_date[6:]))

    # if basetime is 00
    # assign starting date minus 1 day
    # since we need the 12 hours upfront
    # (the day before from 12 UTC to current day 00 UTC)
        start = start - datetime.timedelta(days=1)

    print('Prepare ' + start.strftime("%Y%m%d") +
          '/to/' + end.strftime("%Y%m%d"))

    # create output dir if necessary
    if not os.path.exists(c.outputdir):
        make_dir(c.outputdir)

    # get all files with flux data to be deaccumulated
    inputfiles = UioFiles(c.inputdir, '*OG_acc_SL*.' + str(c.ppid) + '.*')

    # deaccumulate the flux data
    flexpart = EcFlexpart(c, fluxes=True)
    flexpart.write_namelist(c)
    flexpart.deacc_fluxes(inputfiles, c)

    # get a list of all other files
    inputfiles = UioFiles(c.inputdir, '????__??.*' + str(c.ppid) + '.*')

    # produce FLEXPART-ready GRIB files and process them -
    # copy/transfer/interpolate them or make them GRIB2
    flexpart = EcFlexpart(c, fluxes=False)
    flexpart.create(inputfiles, c)
    if c.stream.lower() == 'elda' and c.doubleelda:
        flexpart.calc_extra_elda(c.inputdir, c.prefix)
    flexpart.process_output(c)

    # check if in debugging mode, then store all files
    # otherwise delete temporary files
    if c.debug:
        print('\nTemporary files left intact')
    else:
        clean_up(c)

    return

    main()

# @Author: Anne Philipp (University of Vienna)
#
# @Date: May 2018
#
#    October 2014 - Anne Fouilloux (University of Oslo)
#        - created functions silent_remove and product (taken from ECMWF)
#
#    November 2015 - Leopold Haimberger (University of Vienna)
#        - created functions: interpret_args_and_control, clean_up
#          my_error, normal_exit, init128, to_param_id
#
#    April - December 2018 - Anne Philipp (University of Vienna):
#        - moved all non class methods from former file Flexparttools in here
#        - separated args and control interpretation
#        - added functions get_list_as_string, read_ecenv, send_mail, make_dir,
#          put_file_to_ecserver, submit_job_to_ecserver, get_informations,
#          get_dimensions, execute_subprocess, none_or_int, none_or_str
#
#
#
#
#    none_or_str
#    none_or_int
#    get_cmdline_args
#    read_ecenv
#    clean_up
#    my_error
#    send_mail
#    normal_exit
#    product
#    silent_remove
#    init128
#    to_param_id
#    get_list_as_string
#    make_dir
#    put_file_to_ecserver
#    submit_job_to_ecserver
#    get_informations
#    get_dimensions
#    execute_subprocess
'''This module contains a collection of diverse tasks within flex_extract.
'''

# MODULES

import os
import errno
import sys
# pylint: disable=unused-import
try:
    import exceptions
    import builtins as exceptions
# pylint: enable=unused-import

# METHODS

def setup_controldata():
    '''Collects, stores and checks controlling arguments from command line,
    CONTROL file and ECMWF_ENV file.

    Return
    ------

    ppid : str
        Parent process id.

    queue : str

    job_template : str
        Name of the job template file for submission to ECMWF server.
    '''
    import _config
    from Classes.ControlFile import ControlFile

    args = get_cmdline_args()
    if os.path.isfile(_config.PATH_ECMWF_ENV):
        env_parameter = read_ecenv(_config.PATH_ECMWF_ENV)
        c.assign_envs_to_control(env_parameter)
    c.check_conditions(args.queue)

    return c, args.ppid, args.queue, args.job_template

def none_or_str(value):
    '''Converts the input string into Pythons None type if it
    contains the string "None".

    value : str
        String to be checked for the "None" word.

    Return
    ------
    None or value:
        Return depends on the content of the input value. If it was "None",
        then the Python type None is returned, otherwise the string itself.
    '''
    if value == 'None':
        return None
    return value

def none_or_int(value):
    '''Converts the input string into Pythons None-type if it
    contains string "None"; otherwise it is converted to an integer value.

    value : str

    Return
    ------
    None or int(value):
        then the python type None is returned. Otherwise the string is
        converted into an integer value.
    '''
    return int(value)

def get_cmdline_args():

    Return
    ------
        Contains the command line arguments from the script / program call.
    '''

    parser = ArgumentParser(description='Retrieve FLEXPART input from \
                                ECMWF MARS archive',

    # control parameters that override control file values
    parser.add_argument("--start_date", dest="start_date",
                        help="start date YYYYMMDD")
    parser.add_argument("--end_date", dest="end_date",
                        help="end_date YYYYMMDD")
    parser.add_argument("--date_chunk", dest="date_chunk",
                        type=none_or_int, default=None,
                        help="# of days to be retrieved at once")
    parser.add_argument("--job_chunk", dest="job_chunk",
                        help="# of days to be retrieved within a single job")
                        help="The file with all CONTROL parameters.")
    parser.add_argument("--basetime", dest="basetime",
                        help="base time such as 0 or 12 (for half day retrievals)")
    parser.add_argument("--step", dest="step",
                        help="Forecast steps such as 00/to/48")
    parser.add_argument("--levelist", dest="levelist",
                        help="Vertical levels to be retrieved, e.g. 30/to/60")
    parser.add_argument("--area", dest="area",
                        help="area, defined by north/west/south/east")

    # some switches
    parser.add_argument("--debug", dest="debug",
                        help="debug mode - temporary files will be conserved")
    parser.add_argument("--oper", dest="oper",
                        help='operational mode - prepares dates from '
                        'environment variables')
    parser.add_argument("--request", dest="request",
                        help="list all MARS requests in file mars_requests.dat")
    parser.add_argument("--public", dest="public",
                        help="public mode - retrieves public datasets")
    parser.add_argument("--rrint", dest="rrint",
                        help='Selection of old or new  '
                        'interpolation method for precipitation:\n'
                        '     0 - old method\n'
                        '     1 - new method (additional subgrid points)')

    # set directories
    parser.add_argument("--inputdir", dest="inputdir",
                        help='Path to temporary directory for '
                        'retrieved grib files and other processing files.')
    parser.add_argument("--outputdir", dest="outputdir",
                        help='Path to final directory where '
                        'FLEXPART input files will be stored.')

    # this is only used by prepare_flexpart.py to rerun a postprocessing step
    parser.add_argument("--ppid", dest="ppid",
                        help='This is the specify the parent process id of a '
                        'single flex_extract run to identify the files. '
                        'It is the second number in the GRIB files.')

                        type=none_or_str, default="job.temp",
                        help='Job template file. Will be used for submission '
                        'to the batch system on the ECMWF server after '
                        'modification.')
    parser.add_argument("--queue", dest="queue",
                        help='The name of the ECMWF server name where the'
                        'job script is to be submitted ' 
                        '(e.g. ecgate | cca | ccb)')

def read_ecenv(filepath):
    '''Reads the file into a dictionary where the key values are the parameter
    names.

    filepath : str
        Path to file where the ECMWF environment parameters are stored.

    Return
    ------
    envs : dict
        Contains the environment parameter ecuid, ecgid, gateway
        and destination for ECMWF server environments.
    '''
    envs = {}
    try:
        with open(filepath, 'r') as f:
            for line in f:
                data = line.strip().split()
                envs[str(data[0])] = str(data[1])

        sys.exit('\n... Error occured while trying to read ECMWF_ENV '
                 'file: ' + str(filepath))

    return envs

def clean_up(c):
    '''Remove files from the intermediate directory (inputdir).

    It keeps the final FLEXPART input files if program runs without
    ECMWF API and keywords "ectrans" or "ecstorage" are set to "1".

    Return
    ------

    '''

    print("... clean inputdir!")

    cleanlist = [filename for filename in
                 glob.glob(os.path.join(c.inputdir, "*"))
                 if not os.path.basename(filename).startswith(c.prefix)]

    if cleanlist:
        for element in cleanlist:
            silent_remove(element)
        print("... done!")
    else:
        print("... nothing to clean!")

    return

def my_error(message='ERROR'):
    '''Prints a specified error message which can be passed to the function
    before exiting the program.

    message : str, optional
        Error message. Default value is "ERROR".

    Return
    ------

    '''

    trace = '\n'.join(traceback.format_stack())
    full_message = message + '\n\n' + trace

    print(full_message)

    sys.exit(1)

    return

def send_mail(users, success_mode, message):
    '''Prints a specific exit message which can be passed to the function.

    users : list of str
        Contains all email addresses which should be notified.
        It might also contain just the ecmwf user name which wil trigger
        mailing to the associated email address for this user.

    success_mode : str
        States the exit mode of the program to put into
        the mail subject line.

        Message for exiting program. Default value is "Done!".

    Return
    ------

    '''

    for user in users:
        if '${USER}' in user:
            user = os.getenv('USER')
        try:
            p = subprocess.Popen(['mail', '-s flex_extract_v7.1 ' +
                                  success_mode, os.path.expandvars(user)],
            pout = p.communicate(input=message + '\n\n')[0]
        except ValueError as e:
            print('... ERROR: ' + str(e))
            sys.exit('... Email could not be sent!')
        else:
            print('Email sent to ' + os.path.expandvars(user))

    return

def normal_exit(message='Done!'):

    Return
    ------

    '''

    print(str(message))

    return

def product(*args, **kwds):
    '''Creates combinations of all passed arguments.

    This method combines the single characters of the passed arguments
    with each other in a way that each character of each argument value
    will be combined with each character of the other arguments as a tuple.

    Note
    ----
    This method is taken from an example at the ECMWF wiki website.
    https://software.ecmwf.int/wiki/display/GRIB/index.py; 2018-03-16

    It was released under the following license:
    https://confluence.ecmwf.int/display/ECC/License

    Example
    -------
    product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy

    product(range(2), repeat = 3) --> 000 001 010 011 100 101 110 111

    \*args : list or str
        Positional arguments (arbitrary number).

    \*\*kwds : dict
        Contains all the keyword arguments from \*args.

    Return
    ------
    prod : :obj:`tuple`
        Return will be done with "yield". A tuple of combined arguments.
        See example in description above.
    '''
    try:
        pools = [tuple(arg) for arg in args] * kwds.get('repeat', 1)
        result = [[]]
        for pool in pools:
            result = [x + [y] for x in result for y in pool]
        for prod in result:
            yield tuple(prod)
    except TypeError as e:
        sys.exit('... PRODUCT GENERATION FAILED!')

    return

def silent_remove(filename):
    '''Remove file if it exists.
    The function does not fail if the file does not exist.

    filename : str
        The name of the file to be removed without notification.

    Return
    ------

    '''
    try:
        os.remove(filename)
        # errno.ENOENT  =  no such file or directory
        if e.errno == errno.ENOENT:
            pass
        else:
            raise  # re-raise exception if a different error occured

    return

def init128(filepath):
    '''Opens and reads the grib file with table 128 information.

        Path to file of ECMWF grib table number 128.

    Return
    ------
    table128 : dict
        Contains the ECMWF grib table 128 information.
        The key is the parameter number and the value is the
        short name of the parameter.
    '''
    table128 = dict()
    try:
        with open(filepath) as f:
            fdata = f.read().split('\n')

        sys.exit('\n... Error occured while trying to read parameter '
                 'table file: ' + str(filepath))
    else:
        for data in fdata:
            if data != '' and data[0] != '!':
                table128[data[0:3]] = data[59:65].strip()

    return table128

def to_param_id(pars, table):
    '''Transform parameter names to parameter ids with ECMWF grib table 128.

    pars : str
        Addpar argument from CONTROL file in the format of
        parameter names instead of IDs. The parameter short
        names are separated by "/" and passed as one single string.

    table : dict

    Return
    ------
    ipar : list of int
        List of addpar parameters from CONTROL file transformed to
        parameter ids in the format of integer.
    '''
    if not pars:
        return []
    if not isinstance(pars, str):
        pars = str(pars)

    cpar = pars.upper().split('/')
    ipar = []
    for par in cpar:
        par = par.strip()
        for k, v in table.items():
            if par.isdigit():
                par = str(int(par)).zfill(3)
            if par == k or par == v:
                ipar.append(int(k))
                break
        else:
            print('Warning: par ' + par + ' not found in table 128')

    return ipar

def to_param_id_with_tablenumber(pars, table):
    '''Transform parameter names to parameter IDs and add table ID.

    Conversion with ECMWF grib table 128.

    pars : str
        parameter names instead of ID. The parameter short

    Return
    ------
    spar : str
        parameter IDs in the format of integer.
    '''
        return []

    spar = []
                spar.append(k + '.128')
                break
        else:
            print('\n\n\t\tWarning: par ' + par + ' not found in table 128\n\n')

    return '/'.join(spar)

def get_list_as_string(list_obj, concatenate_sign=', '):
    '''Converts a list of arbitrary content into a single string using a given
    concatenation character.

    list_obj : list of *
        A list with arbitrary content.

    concatenate_sign : str, optional
        A string which is used to concatenate the single
        list elements. Default value is ", ".

    Return
    ------
    str_of_list : str
        The content of the list as a single string.
    '''

    if not isinstance(list_obj, list):
        list_obj = list(list_obj)
    str_of_list = concatenate_sign.join(str(l) for l in list_obj)

    return str_of_list

def make_dir(directory):
    '''Creates a directory.

    If the directory already exists, an information is printed and the creation 
    skipped. The program stops only if there is another problem.

    directory : str
        The path to directory which should be created.

    Return
    ------

    '''
    try:
        os.makedirs(directory)
        # errno.EEXIST = directory already exists
        if e.errno == errno.EEXIST:
            print('INFORMATION: Directory {0} already exists!'.format(directory))
        else:
            raise # re-raise exception if a different error occured

    return

def put_file_to_ecserver(ecd, filename, target, ecuid, ecgid):
    '''Uses the ecaccess-file-put command to send a file to the ECMWF servers.

    Note
    ----
    The return value is just for testing reasons. It does not have
    to be used from the calling function since the whole error handling
    is done in here.

    ecd : str
        The path were the file is stored.

        The name of the file to send to the ECMWF server.

        The target queue where the file should be sent to.

    ecuid : str

    ecgid : str

    Return
    ------

    '''

    try:
        subprocess.check_output(['ecaccess-file-put',
                                 ecd + '/' + filename,
                                 target + ':/home/ms/' +
                                 ecgid + '/' + ecuid +
                                 '/' + filename],
                                stderr=subprocess.STDOUT)
    except subprocess.CalledProcessError as e:
        print('... ERROR CODE: ' + str(e.returncode))
        print('... ERROR MESSAGE:\n \t ' + str(e))

        print('\n... Do you have a valid ecaccess certification key?')
        sys.exit('... ECACCESS-FILE-PUT FAILED!')

        print('\n... Most likely the ECACCESS library is not available!')

    return

def submit_job_to_ecserver(target, jobname):
    '''Uses ecaccess-job-submit command to submit a job to the ECMWF server.

    Note
    ----

        The target where the file should be sent to, e.g. the queue.

    jobname : str
        The name of the jobfile to be submitted to the ECMWF server.

    Return
    ------
    job_id : int
        The id number of the job as a reference at the ECMWF server.
    '''

    try:
        job_id = subprocess.check_output(['ecaccess-job-submit', '-queueName',
                                          target, jobname])

        sys.exit('... ecaccess-job-submit FAILED!')

    return job_id.decode()

def get_informations(filename):
    '''Extracts basic information from a sample grib file.

    This information is needed for later use and the
    initialization of numpy arrays where data are stored.

            Name of the file which will be opened to extract basic information.

    Return
    ------
    data : dict
        Contains basic informations of the ECMWF grib files, e.g.
        'Ni', 'Nj', 'latitudeOfFirstGridPointInDegrees',
        'longitudeOfFirstGridPointInDegrees', 'latitudeOfLastGridPointInDegrees',
        'longitudeOfLastGridPointInDegrees', 'jDirectionIncrementInDegrees',
        'iDirectionIncrementInDegrees', 'missingValue'
    '''
    from eccodes import codes_grib_new_from_file, codes_get, codes_release

    data = {}

    # --- open file ---
    print("Opening grib file for extraction of information --- %s" % filename)
    with open(filename, 'rb') as f:
        # load first message from file
        gid = codes_grib_new_from_file(f)

        # information needed from grib message
        keys = ['Ni',
                'Nj',
                'latitudeOfFirstGridPointInDegrees',
                'longitudeOfFirstGridPointInDegrees',
                'latitudeOfLastGridPointInDegrees',
                'longitudeOfLastGridPointInDegrees',
                'jDirectionIncrementInDegrees',
                'iDirectionIncrementInDegrees',
                'missingValue',
               ]

        print('\nInformation extracted: ')
        for key in keys:
            # Get the value of the key in a grib message.
            data[key] = codes_get(gid, key)
            print("%s = %s" % (key, data[key]))

        # Free the memory for the message referred as gribid.
        codes_release(gid)

    return data

def get_dimensions(info, purefc, dtime, index_vals, start_date, end_date):
    '''This function specifies the correct dimensions for x, y, and t.

    info : dict

    dtime : str

    index_vals : list of list of str
        Contains the values from the keys used for a distinct selection
        of GRIB messages in processing the grib files.
        Content looks like e.g.:
        index_vals[0]: ('20171106', '20171107', '20171108') ; date
        index_vals[1]: ('0', '1200', '1800', '600') ; time
        index_vals[2]: ('0', '12', '3', '6', '9') ; stepRange

        The start date of the retrieval job.

    end_date : str
        The end date of the retrieval job.

    Return
    ------
    (ix, jy, it) : tuple of int
        Dimension in x-direction, y-direction and in time.
    '''

    ix = info['Ni']

    jy = info['Nj']

    if not purefc:
        it = ((end_date - start_date).days + 1) * 24 // int(dtime)
    else:
        # #no of step * #no of times * #no of days
        it = len(index_vals[2]) * len(index_vals[1]) * len(index_vals[0])

    return (ix, jy, it)

def execute_subprocess(cmd_list, error_msg='SUBPROCESS FAILED!'):
    '''Executes a command via a subprocess.

    Error handling is done if an error occures.

    cmd_list : list of str
        A list of the components for the command line execution. 
        They will be concatenated with blank space for the command 
        to be submitted, like ['mv', file1, file2] for mv file1 file2.

    Return
    ------
    error_msg : str, optional
        Error message if the subprocess fails.
        By default it will just say "SUBPROCESS FAILED!".
    '''

    try:
        subprocess.check_call(cmd_list)

        sys.exit('... ' + error_msg)

    return

def generate_retrieval_period_boundary(c):
    '''Generates retrieval period boundary datetimes from CONTROL information.

    Return
    ------
    start_period : datetime
        The first timestamp of the actual retrieval period disregarding
        the temporary times which were used for processing reasons.

    end_period : datetime
        The last timestamp of the actual retrieval period disregarding
    '''
    # generate start and end timestamp of the retrieval period
    start_period = datetime.strptime(c.start_date + c.time[0], '%Y%m%d%H')
    start_period = start_period + timedelta(hours=int(c.step[0]))
    end_period = datetime.strptime(c.end_date + c.time[-1], '%Y%m%d%H')
    end_period = end_period + timedelta(hours=int(c.step[-1]))

    return start_period, end_period

#
# @Date: November 2018
#
#
#
#
'''This module contains check methods for the CONTROL paramaters.
'''

# MODULES

import os
import sys
from datetime import datetime
try:

from Mods.tools import my_error, silent_remove

def check_logicals_type(c, logicals):
    '''Check that the logical variables have correct type integer.

    logicals : list of (str or int)
        Names of the switches that are used to control the flow of the
        program.

    Return
    ------

    '''

    for var in logicals:
        if not isinstance(getattr(c, var), int):
            setattr(c, var, int(getattr(c, var)))

    return

def check_grid(grid):
    '''Convert grid into correct Lat/Lon format. E.g. '0.5/0.5'

    Checks on format of original grid. Wether it is in the order of 1000 or 1.
    Convert to correct grid format and substitute into "Lat/Lon" format string.

    grid : str
        Contains grid information

    Return
    ------
    grid : str
        Contains grid in format Lat/lon. E.g. 0.1/0.1
    '''

    if 'N' in grid:
        return grid
    if '/' in grid:
        gridx, gridy = grid.split('/')
        if gridx == gridy:
            grid = gridx
        else:
            raise ValueError('GRID parameter contains two '
                             'different values: %s' (grid))
    # # determine grid format
    # if float(grid) / 100. >= 0.5:
    #    # grid is defined in 1/1000 degrees; old format
    #    grid = '{}/{}'.format(float(grid) / 1000.,
    #                          float(grid) / 1000.)
    # elif float(grid) / 100. < 0.5:
    #    # grid is defined in normal degree; new format
    #    grid = '{}/{}'.format(float(grid), float(grid))

    # determine grid format
    # assumes that nobody wants grid spacings of 20 deg or more
    if float(grid) >= 20.:
        # grid is defined in 1/1000 degree; old format
        grid = '{}/{}'.format(float(grid) / 1000., float(grid) / 1000.)
    else:
        # grid is defined in degree; new format
        grid = '{}/{}'.format(float(grid), float(grid))

    return grid

def check_area(grid, area, upper, lower, left, right):
    '''Defines the correct area string.

    Checks the format of the four area components wether it is on
    the order of 1000 or 1. 
    Also checks wether area was already set on command line, 
    then the four components are overwritten.
    Convert to correct format of the order of magnitude "1" and sets the
    area parameter (North/West/South/East).
    E.g.: -5./20./10./10.

    grid : str
        Contains grid information.

    area : str
        Contains area informtion.

    upper : str
        The northernmost latitude.

    lower : str
        The southernmost latitude.

    left : str
        The westernmost longitude.

    right : str
        The easternmost longitude.

    Return
    ------
    grid : str
        Contains grid in format lat/lon. E.g. 0.1/0.1
    '''
    if 'N' in grid:  # Gaussian output grid
        area = 'G'
        return area

    # if area was provided decompose area into its 4 components
    if area:
        components = area.split('/')
        upper, left, lower, right = components

    # determine area format
    if all([(abs(float(upper) / 10000.) >= 0.01 or float(upper) / 1000. == 0.),
            (abs(float(lower) / 10000.) >= 0.01 or float(lower) / 1000. == 0.),
            (abs(float(left) / 10000.) >= 0.01 or float(left) / 1000. == 0.),
            (abs(float(right) / 10000.) >= 0.01 or float(right) / 1000. == 0.)]):
        # area is defined in 1/1000 degrees; old format
        area = '{}/{}/{}/{}'.format(float(upper) / 1000.,
                                    float(left) / 1000.,
                                    float(lower) / 1000.,
                                    float(right) / 1000.)
    elif all([abs(float(upper) / 10000.) < 0.05,
              abs(float(lower) / 10000.) < 0.05,
              abs(float(left) / 10000.) < 0.05,
              abs(float(right) / 10000.) < 0.05]):
        # area is already in new format
        area = '{}/{}/{}/{}'.format(float(upper),
                                    float(left),
                                    float(lower),
                                    float(right))
    else:
        raise ValueError('Area components have inconsisten or unrecognised '
                         'formats (upper, lower, left, right): '
                         '{}/{}/{}/{}'.format(str(upper), str(lower),
                                              str(left), str(right)))

    return area

def check_levels(levelist, level):
    '''Defines correct level list and guarantees that the maximum level is
    one of the available maximum levels.

        Specifies the level list.
        Examples: model level: 1/to/137, pressure levels: 500/to/1000

    level : str

    Return
    ------
        Specifies the required levels. It has to have a valid
        correspondence to the selected levtype.

    level : str
        Specifies the maximum level. It has to be one of the
        available maximum level numbers as contained in the variable
        MAX_LEVEL_LIST in "_config": [16, 19, 31, 40, 50, 60, 62, 91, 137]

    '''
    # assure consistency of levelist and level
    if not levelist and not level:
        raise ValueError('ERROR: neither levelist nor level '
                         'specified in CONTROL file')
    elif not levelist and level:
        levelist = '1/to/' + level
    elif (levelist and not level) or \
         (levelist[-1] != level[-1]):
        level = levelist.split('/')[-1]
    else:
        pass

    # check if max level is a valid level
    if int(level) not in _config.MAX_LEVEL_LIST:
        raise ValueError('ERROR: \n'
                         'LEVEL must be the maximum level of a specified '
                         'level list from ECMWF, e.g. {} \n'
                         'Check parameter "LEVEL" or the max level of '
                         '"LEVELIST"!'.format(str(_config.MAX_LEVEL_LIST)))

    return levelist, level

def check_ppid(c, ppid):
    '''Sets the current PPID.

    ppid : int or None
        Contains the ppid number provided by the command line parameter
        of is None otherwise.

    Return
    ------

    '''

    if not ppid:
        c.ppid = str(os.getppid())
    else:
        c.ppid = ppid

    return

def check_purefc(ftype):
    '''Check for a pure forecast mode.

    ftype : list of str

    Return
    ------
    True or False:
        True if pure forecasts are to be retrieved. False if there are
        analysis fields in between.
    '''

    if 'AN' not in ftype and '4V' not in ftype:
        # pure forecast
        return 1

    return 0

def check_step(step):
    '''Checks on step format and convert into a list of steps.

    If the steps were defined with "to" and "by" they are converted into
    a list of steps. If the steps were set in a string, it is
    converted into a list.

        Valid values are hours (HH) from forecast base time.

    Return
    ------
    step : list of str
        List of forecast steps in format e.g. [001, 002, ...]
    '''
    import numpy as np

    if '/' in step:
        steps = step.split('/')
        if 'to' in step.lower() and 'by' in step.lower():
            ilist = np.arange(int(steps[0]),
                              int(steps[2]) + 1,
                              int(steps[4]))
            step = ['{:0>3}'.format(i) for i in ilist]
        elif 'to' in step.lower() and 'by' not in step.lower():
            my_error(step + ':\n' +
                     'if "to" is used in steps parameter, '
                     'please use "by" as well')
        else:
            step = steps

    if not isinstance(step, list):
        step = [step]

    return step

def check_type(ftype, steps):
    '''Check if type variable is of type list and if analysis field has
    forecast step 0.

    ftype : list of str or str

    steps : str

    Return
    ------
    '''
    if not isinstance(ftype, list):
        ftype = [ftype]

    for i, val in enumerate(ftype):
        if ftype[i] == 'AN' and int(steps[i]) != 0:
            print('For analysis data, STEP = 0 is needed. Setting to 0 now.)')
            ftype[i] = 0

    return ftype

def check_time(ftime):
    '''Check if time variable is of type list. Otherwise convert to list.

    ftime : list of str or str

    Return
    ------
    ftime : list of str
    '''
    if not isinstance(ftime, list):
        ftime = [ftime]

    return ftime

def check_len_type_time_step(ftype, ftime, steps, maxstep, purefc):
    '''Check if

    steps : str

        The maximum forecast time step in hours from the forecast base time.
        This is the maximum step for non-flux (not accumulated) forecast data.

    Return
    ------

    steps : str
    '''
    if not len(ftype) == len(ftime) == len(steps):
        raise ValueError('ERROR: The number of field types, times and steps '
                         'are not the same! Please check the settings in the '
                         'CONTROL file!')

    # if pure forecast is selected and only one field type/time is set
    # prepare a complete list of type/time/step combination upto maxstep
    if len(ftype) == 1 and purefc:
        nftype = []
        nsteps = []
        nftime = []
        for i in range(0, maxstep + 1):
            nftype.append(ftype[0])
            nsteps.append('{:0>3}'.format(i))
            nftime.append(ftime[0])
        return nftype, nftime, nsteps

    return ftype, ftime, steps

def check_mail(mail):
    '''Check the string of mail addresses, separate them and convert to list.

    mail : list of str or str
        Contains email addresses for notifications.
        It might also contain just the ecmwf user name which will trigger

    Return
    ------
    mail : list of str

    '''
    if not isinstance(mail, list):
        if ',' in mail:
            mail = mail.split(',')
        elif ' ' in mail:
            mail = mail.split()
        else:
            mail = [mail]

    return mail

def check_queue(queue, gateway, destination, ecuid, ecgid):
    '''Check if the necessary ECMWF parameters are set if the queue is
    one of the QUEUES_LIST (in _config).

    queue : str
        Name of the queue if submitted to the ECMWF servers.
        Used to check if ecuid, ecgid, gateway and destination
        are set correctly and are not empty.

        The address of the gateway server.

        The name of the destination of the gateway server for data
        transfer through ectrans. E.g. name@genericSftp

    ecuid : str
        ECMWF user id.

    ecgid : str
        ECMWF group id.

    Return
    ------

    '''
    if queue in _config.QUEUES_LIST and \
            (not ecuid or not ecgid):
        raise ValueError('\nEnvironment variables ECUID '
                         'and ECGID were not set properly! \n '
                         'Please check for existence of file "ECMWF_ENV" '
                         'in the run directory!')
    elif queue in _config.QUEUES_LIST and \
             (not gateway or not destination):
        print('WARNING: Parameters GATEWAY and DESTINATION were '
              'not properly set for working on ECMWF server. \n'
              'There will be no transfer of output files to the '
              'local gateway server possible!')
    return

def check_pathes(idir, odir, fpdir, fedir):
    '''Check if output and flexpart pathes are set.

    idir : str
        Path to the temporary directory for MARS retrieval data.

    odir : str
        Path to the final output directory where the FLEXPART input files
        will be stored.

    fpdir : str
        Path to FLEXPART root directory.

    fedir : str
        Path to flex_extract root directory.

    Return
    ------
    odir : str

    fpdir : str

    '''
    if not fpdir:
        fpdir = fedir

    if not odir:
        odir = idir

    return odir, fpdir

def check_dates(start, end):
    '''Checks if there is at least a start date for a one-day retrieval.

    Checks if end date lies after start date and end date is set.

    start : str

    end : str

    Return
    ------
    start : str

    end : str

    '''
    # check for having at least a starting date
    # otherwise program is not allowed to run
    if not start:
        raise ValueError('start_date was neither specified in command line nor '
                         'in CONTROL file.\n'
                         'Try "{} -h" to print usage information'
                         .format(sys.argv[0].split('/')[-1]))

    # retrieve just one day if end_date isn't set
    if not end:
        end = start

    dstart = datetime.strptime(start, '%Y%m%d')
    dend = datetime.strptime(end, '%Y%m%d')
    if dstart > dend:
        raise ValueError('ERROR: Start date is after end date! \n'
                         'Please adapt the dates in CONTROL file or '
                         'command line! (start={}; end={})'.format(start, end))

    return start, end

def check_maxstep(maxstep, steps):
    '''Convert maxstep into integer if it is already given. Otherwise, select
    maxstep by going through the steps list.

    maxstep : str
        This is the maximum step for non flux (accumulated) forecast data.

    steps : str

    Return
    ------

    '''
    # if maxstep wasn't provided
    # search for it in the "step" parameter
    if not maxstep:
        maxstep = 0
        for s in steps:
            if int(s) > maxstep:
                maxstep = int(s)
    else:
        maxstep = int(maxstep)

    return maxstep

def check_basetime(basetime):
    '''Check if basetime is set and contains one of the two
    possible values (0, 12).

    basetime : int or str or None

    Return
    ------
    basetime : int or None
    '''
    if basetime is not None:
        basetime = int(basetime)
        if basetime != 0 and basetime != 12:
            raise ValueError('ERROR: Basetime has an invalid value '
                             '-> {}'.format(str(basetime)))
    return basetime

def check_request(request, marsfile):
    '''Check if there is an old MARS request file; if so, remove it.

        Selects the mode of retrieval.
        0: Retrieves the data from ECMWF.
        1: Prints the mars requests to an output file.
        2: Retrieves the data and prints the mars request.

        Path to the mars request file.

    Return
    ------

    '''
    if request != 0:
        if os.path.isfile(marsfile):
            silent_remove(marsfile)
    return

def check_public(public, dataset):
    '''Check wether the dataset parameter is set to a
    public data set.

        Specifies if public data are to be retrieved or not.

        Specific name which identifies the public dataset.

    Return
    ------

    '''
    if public and not dataset:
        raise ValueError('ERROR: If public MARS data are to be retrieved, '
                         'the "dataset"-parameter has to be set, too!')
    return

def check_acctype(acctype, ftype):
    '''Guarantees that the accumulation field type is set.

    If not set, it is derived with the old method (TYPE[1]).

    Return
    ------
    '''
    if not acctype:
        print('... Control parameter ACCTYPE was not defined.')
        try:
            if len(ftype) == 1 and ftype[0] != 'AN':
                print('... Use same field type as for the non-flux fields.')
                acctype = ftype[0]
            elif len(ftype) > 1 and ftype[1] != 'AN':
                print('... Use old setting by using TYPE[1] for flux forecast!')
                acctype = ftype[1]
        except:
            raise ValueError('ERROR: Accumulation field type could not be set!')
    else:
        if acctype.upper() == 'AN':
            raise ValueError('ERROR: Accumulation forecast fields can not be '
                             'of type "analysis"!')
    return acctype

def check_acctime(acctime, marsclass, purefc, time):
    '''Guarantees that the accumulation forecast times were set.

    If not set, setting the value to some of the most commonly used data sets
    is attempted. Otherwise, an eror is raised.

        The starting time for the accumulated forecasts.

        ECMWF data classification identifier.

    Return
    ------
    '''

    if not acctime:
        print('... Control parameter ACCTIME was not set.')
        print('... Value will be set depending on field type:\n '
              '\t\t EA=06/18\n\t\t EI/OD=00/12\n\t\t EP=18')
        if marsclass.upper() == 'EA': # Era 5
            acctime = '06/18'
        elif marsclass.upper() == 'EI': # Era-Interim
            acctime = '00/12'
        elif marsclass.upper() == 'EP': # CERA
            acctime = '18'
        elif marsclass.upper() == 'OD' and not purefc: # On-demand
        elif marsclass.upper() == 'OD' and purefc: # On-demand
            acctime = time[0]
        else:
            raise ValueError('ERROR: Accumulation forecast time can not be'
                             'derived automatically!')
    return acctime

def check_accmaxstep(accmaxstep, marsclass, purefc, maxstep):
    '''Guarantees that the accumulation forecast step was set.

        This is the maximum step for non-flux (accumulated) forecast data.

    Return
    ------
    '''
    if not accmaxstep:
        print('... Control parameter ACCMAXSTEP was not set.')
        print('... Value will be set depending on field type/time: '
              '\n\t\t EA/EI/OD=12\n\t\t EP=24')
        if marsclass.upper() in ['EA', 'EI', 'OD'] and not purefc:
            # Era 5, Era-Interim, On-demand operational
            accmaxstep = '12'
            accmaxstep = '24'
        elif purefc and accmaxstep != maxstep:
            accmaxstep = maxstep
            print('... For pure forecast mode, the accumulated forecast must '
                  'have the same maxstep as the normal forecast fields!\n'
                  '\t\t ACCMAXSTEP was set to MAXSTEP!')
        else:
            raise ValueError('ERROR: Accumulation forecast step can not be'
    else:
        if purefc and int(accmaxstep) != int(maxstep):
    return accmaxstep

def check_addpar(addpar):
    '''Check that addpar has correct format of additional parameters in
    a single string, so that it can be easily appended to the hard-coded
    parameters retrieved in any case.

    addpar : str or list of str

    Return
    ------
    '''

    if addpar and isinstance(addpar, str):
        if '/' in addpar:
            parlist = addpar.split('/')
            parlist = [p for p in parlist if p is not '']
        else:
            parlist = [addpar]

        addpar = '/' + '/'.join(parlist)

    return addpar

def check_job_chunk(job_chunk):
    '''Checks that if job chunk is set, the number is positive and nonzero.

    job_chunk : int
        The number of days for a single job script.

    Return
    ------
    '''
    if not job_chunk:
        return job_chunk
    else:
        job_chunk = int(job_chunk)

    if job_chunk < 0:
        raise ValueError('ERROR: The number of job chunks is negative!\n'
                         'It has to be a positive number!')
    elif job_chunk == 0:
        job_chunk = None
    else:
        pass

    return job_chunk

def check_number(number):
    '''Check for correct string format of ensemble member numbers.

        List of ensemble member forecast runs.

    Return
    ------
        String with list of ensemble member forecast runs. E.g. '01/02/03/04'
    '''

    if '/' in number:
        numbers = number.split('/')
        if 'to' in number.lower() and 'by' in number.lower():
            number = '{:0>3}'.format(int(numbers[0])) + '/TO/' + \
                     '{:0>3}'.format(int(numbers[2])) + '/BY/' + \
                     '{:0>3}'.format(int(numbers[4]))
        elif 'to' in number.lower() and 'by' not in number.lower():
                     '{:0>3}'.format(int(numbers[2]))
        else:
            numbers = ['{:0>3}'.format(i) for i in numbers]
            number = '{:0>3}/'.join(numbers)
    elif number.isdigit():
        number = '{:0>3}'.format(int(number))
    else:
        pass

    return number

#
# @Date: March 2018
#
#
#        - migration of the methods dapoly and darain from Fortran
#          (flex_extract_v6 and earlier) to Python
#
#    April 2018 - Anne Philipp (University of Vienna):
#        - outsourced the disaggregation functions dapoly and darain
#          to a new module named disaggregation
#        - added the new disaggregation method for precipitation
#
#
#
#
#    - dapoly
#    - darain
#    - IA3
'''Disaggregation of deaccumulated flux data from an ECMWF model FG field.

Initially the flux data to be concerned are:
    - large-scale precipitation
    - convective precipitation
    - surface sensible heat flux
    - surface solar radiation
    - u stress
    - v stress

Different versions of disaggregation is provided for rainfall
data (darain, modified linear) and the surface fluxes and
stress data (dapoly, cubic polynomial).
'''

# MODULES

def dapoly(alist):
    """Cubic polynomial interpolation of deaccumulated fluxes.

    Interpolation of deaccumulated fluxes of an ECMWF model FG field
    using a cubic polynomial solution which conserves the integrals
    of the fluxes within each timespan.
    Disaggregation is done for 4 accumluated timespans which
    generates a new, disaggregated value which is output at the
    central point of the 4 accumulation timespans.
    This new point is used for linear interpolation of the complete
    timeseries afterwards.

    alist : list of array of float
        List of 4 timespans as 2-dimensional, horizontal fields.
        E.g. [[array_t1], [array_t2], [array_t3], [array_t4]]

    Return
    ------
    nfield : array of float
        Interpolated flux at central point of accumulation timespan.

    Note
    ----
    March 2000    : P. JAMES
        Original author

    June 2003     : A. BECK
        Adaptations

    November 2015 : Leopold Haimberger (University of Vienna)
        Migration from Fortran to Python

    """

    pya = (alist[3] - alist[0] + 3. * (alist[1] - alist[2])) / 6.
    pyb = (alist[2] + alist[0]) / 2. - alist[1] - 9. * pya / 2.
    pyc = alist[1] - alist[0] - 7. * pya / 2. - 2. * pyb
    pyd = alist[0] - pya / 4. - pyb / 3. - pyc / 2.
    nfield = 8. * pya + 4. * pyb + 2. * pyc + pyd

    return nfield

def darain(alist):
    """Linear interpolation of deaccumulated fluxes.

    Interpolation of deaccumulated fluxes of an ECMWF model FG rainfall
    field using a modified linear solution which conserves the integrals
    Disaggregation is done for 4 accumluated timespans which generates
    a new, disaggregated value which is output at the central point
    of the 4 accumulation timespans. This new point is used for linear
    interpolation of the complete timeseries afterwards.

    Return
    ------

    Note
    ----

    """

    xa = alist[0]
    xb = alist[1]
    xc = alist[2]
    xd = alist[3]
    xa[xa < 0.] = 0.
    xb[xb < 0.] = 0.
    xc[xc < 0.] = 0.
    xd[xd < 0.] = 0.

    xac = 0.5 * xb
    mask = xa + xc > 0.
    xac[mask] = xb[mask] * xc[mask] / (xa[mask] + xc[mask])
    xbd = 0.5 * xc
    mask = xb + xd > 0.
    xbd[mask] = xb[mask] * xc[mask] / (xb[mask] + xd[mask])
    nfield = xac + xbd

def IA3(g):
    """ Interpolation with a non-negative geometric mean based algorithm.

    The original grid is reconstructed by adding two sampling points in each
    data series interval. This subgrid is used to keep all information during
    the interpolation within the associated interval. Additionally, an advanced
    monotonicity filter is applied to improve the monotonicity properties of
    the series.

    Note
    ----
    (C) Copyright 2017-2019
    Sabine Hittmeir, Anne Philipp, Petra Seibert

    This work is licensed under the Creative Commons Attribution 4.0
    International License. To view a copy of this license, visit
    http://creativecommons.org/licenses/by/4.0/ or send a letter to
    Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

    g : list of float
        Complete data series that will be interpolated having
        the dimension of the original raw series.

    Return
    ------
    f : list of float
        The interpolated data series with additional subgrid points.
        Its dimension is equal to the length of the input data series
        times three.

    References
    For more information see article:
    Hittmeir, S.; Philipp, A.; Seibert, P. (2017): A conservative
    interpolation scheme for extensive quantities with application to the
    Lagrangian particle dispersion model FLEXPART.,
    Geoscientific Model Development
    """

    #######################  variable description #############################
    #                                                                         #
    # i      - index variable for looping over the data series                #
    # g      - input data series                                              #
    # f      - interpolated and filtered data series with additional          #
    #          grid points                                                    #
    # fi     - function value at position i, f_i                              #
    # fi1    - first  sub-grid function value f_i^1                           #
    # fi2    - second sub-grid function value f_i^2                           #
    # fip1   - next function value at position i+1, f_(i+1)                   #
    # dt     - time step                                                      #
    # fmon   - monotonicity filter                                            #
    #                                                                         #
    ###########################################################################

    # time step
    dt = 1.0

    ############### Non-negative Geometric Mean Based Algorithm ###############

    # for the left boundary the following boundary condition is valid:
    # the value at t=0 of the interpolation algorithm coincides with the
    # first data value according to the persistence hypothesis
    f = [g[0]]

    # compute two first sub-grid intervals without monotonicity check
    # go through the data series and extend each interval by two sub-grid
    # points and interpolate the corresponding data values
    # except for the last interval due to boundary conditions
    for i in range(0, 2):

        # as a requirement:
        # if there is a zero data value such that g[i]=0, then the whole
        # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0
        # according to Eq. (6)
        if g[i] == 0.:
            f.extend([0., 0., 0.])

        # otherwise the sub-grid values are calculated and added to the list
        else:
            # temporal save of last value in interpolated list
            # since it is the left boundary and hence the new (fi) value
            fi = f[-1]

            # the value at the end of the interval (fip1) is prescribed by the
            # geometric mean, restricted such that non-negativity is guaranteed
            # according to Eq. (25)
            fip1 = min(3. * g[i], 3. * g[i + 1], np.sqrt(g[i + 1] * g[i]))

            # the function value at the first sub-grid point (fi1) is determined
            # according to the equal area condition with Eq. (19)
            fi1 = 3./2.*g[i]-5./12.*fip1-1./12.*fi

            # the function value at the second sub-grid point (fi2) is determined
            # according Eq. (18)
            fi2 = fi1+1./3.*(fip1-fi)

            # add next interval of interpolated (sub-)grid values
            f.append(fi1)
            f.append(fi2)
            f.append(fip1)

    # compute rest of the data series intervals
    for i in range(2, len(g)-1):

            # apply monotonicity filter for interval before
            # check if there is "M" or "W" shape
            if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5]) == -1 \
               and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4]) == -1 \
               and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3]) == -1:

                # the monotonicity filter corrects the value at (fim1) by
                # substituting (fim1) with (fmon), see Eq. (27), (28) and (29)
                fmon = min(3. * g[i - 2],
                           3. * g[i - 1],
                           np.sqrt(max(0, (18. / 13. * g[i - 2] - 5. / 13. * f[-7]) *
                                       (18. / 13. * g[i - 1] - 5. / 13. * f[-1]))))

                # recomputation of the sub-grid interval values while the
                # interval boundaries (fi) and (fip2) remains unchanged
                # see Eq. (18) and (19)
                f[-4] = fmon
                f[-6] = 3./2.*g[i-2]-5./12.*fmon-1./12.*f[-7]
                f[-5] = f[-6]+(fmon-f[-7])/3.
                f[-3] = 3./2.*g[i-1]-5./12.*f[-1]-1./12.*fmon
                f[-2] = f[-3]+(f[-1]-fmon)/3.

        else:
            fi = f[-1]

                # substituting (fim1) with fmon, see Eq. (27), (28) and (29)

    # separate treatment of the final interval

    # as a requirement:
    # if there is a zero data value such that g[i]=0, then the whole
    # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0
    # according to Eq. (6)
    if g[-1] == 0.:
        # apply monotonicity filter for interval before
        # check if there is "M" or "W" shape
        if     np.sign(f[-5]-f[-6]) * np.sign(f[-4]-f[-5]) == -1 \
           and np.sign(f[-4]-f[-5]) * np.sign(f[-3]-f[-4]) == -1 \
           and np.sign(f[-3]-f[-4]) * np.sign(f[-2]-f[-3]) == -1:

            # the monotonicity filter corrects the value at (fim1) by
            # substituting (fim1) with (fmon), see Eq. (27), (28) and (29)
            fmon = min(3. * g[-3],
                       3. * g[-2],
                       np.sqrt(max(0, (18. / 13. * g[-3] - 5. / 13. * f[-7]) *
                                   (18. / 13. * g[-2] - 5. / 13. * f[-1]))))

            # recomputation of the sub-grid interval values while the
            # interval boundaries (fi) and (fip2) remains unchanged
            # see Eq. (18) and (19)
            f[-4] = fmon
            f[-6] = 3./2.*g[-3]-5./12.*fmon-1./12.*f[-7]
            f[-5] = f[-6]+(fmon-f[-7])/3.
            f[-3] = 3./2.*g[-2]-5./12.*f[-1]-1./12.*fmon
            f[-2] = f[-3]+(f[-1]-fmon)/3.

        f.extend([0., 0., 0.])

    # otherwise the sub-grid values are calculated and added to the list
    # using the persistence hypothesis as boundary condition
    else:
        # temporal save of last value in interpolated list
        # since it is the left boundary and hence the new (fi) value
        fi = f[-1]
        # since last interval in series, last value is also fip1
        fip1 = g[-1]
        # the function value at the first sub-grid point (fi1) is determined
        # according to the equal area condition with Eq. (19)
        fi1 = 3./2.*g[-1]-5./12.*fip1-1./12.*fi
        # the function value at the second sub-grid point (fi2) is determined
        # according Eq. (18)
        fi2 = fi1+dt/3.*(fip1-fi)

                       3. * g[-2],

        # add next interval of interpolated (sub-)grid values
        f.append(fi1)
        f.append(fi2)
        f.append(fip1)

    return f