#!/usr/bin/env python
# -*- coding: utf-8 -*-
#*******************************************************************************
# @Author: Anne Fouilloux (University of Oslo)
#
# @Date: October 2014
#
# @Change History:
#
#   November 2015 - Leopold Haimberger (University of Vienna):
#        - optimized display_info
#        - optimized data_retrieve and seperate between python and shell
#          script call
#
#   February 2018 - Anne Philipp (University of Vienna):
#        - applied PEP8 style guide
#        - added documentation
#        - applied some minor modifications in programming style/structure
#        - added writing of mars request attributes to a csv file
#
# @License:
#    (C) Copyright 2014-2019.
#    Anne Philipp, Leopold Haimberger
#
#    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.

# ------------------------------------------------------------------------------
# MODULES
import os
import sys
import subprocess
import traceback

# software specific classes and modules from flex_extract
sys.path.append('../')
import _config
try:
    ec_api = True
    import ecmwfapi
except ImportError:
    ec_api = False

try:
    cds_api = True
    import cdsapi
    cds_api = False
# 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

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

    public : int
        Decides which Web API Server version is used.

    marsclass : str, optional
        Characterisation of dataset.

    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.

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

    levtype : str, optional
        Denotes type of level.

    levelist : str, optional
        Specifies the required levels.

    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
        Specifies the desired triangular truncation of retrieved data,
        before carrying out any other selected post-processing.

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

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

    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
        The version of the dataset.

    number : str, optional
        Selects the member in ensemble forecast run.

    accuracy : str, optional
        Specifies the number of bits per value to be used in the
        generated GRIB coded fields.

    grid : str, optional
        Specifies the output grid which can be either a Gaussian grid
        or a Latitude/Longitude grid.

    gaussian : str, optional
        This parameter is deprecated and should no longer be used.
        Specifies the desired type of Gaussian grid for the output.

    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

        Parameters
        ----------
        server : ECMWFService or ECMWFDataServer
            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
            Specifies the forecast time step from forecast base time.
            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.server = server
        self.public = public
        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
        attrs = vars(self).copy()

        # iterate through all attributes and print them
        # with their corresponding values
        for item in attrs.items():
            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.iterkeys())))
            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.iteritems():
            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]

        # 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)
                tb = sys.exc_info()[2]
                print(traceback.format_exc())
                sys.exit()

        # MARS request via call in shell
        else:
            request_str = 'ret'
            for key, value in attrs.iteritems():
                request_str = request_str + ',' + key + '=' + str(value)
            request_str += ',target="' + target + '"'
            p = subprocess.Popen(['mars', '-p'],
                                 stdin=subprocess.PIPE,
                                 stdout=subprocess.PIPE,
                                 stderr=subprocess.PIPE,
                                 bufsize=1)
            pout = p.communicate(input=request_str)[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

# @Author: Anne Philipp (University of Vienna)
#
# @Date: March 2018
#
#
#    November 2015 - Leopold Haimberger (University of Vienna):
#        - migration of the methods dapoly and darain from Fortran
#          (flex_extract_v6 and earlier) to Python
#
#    April 2018 - Anne Philipp (University of Vienna):
#        - added structured documentation
#        - outsourced the disaggregation functions dapoly and darain
#          to a new module named disaggregation
#        - added the new disaggregation method for precipitation
#
#
#
# @Methods:
#    - 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

# FUNCTIONS
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.

    Parameters
    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                                            #
    #                                                                         #
    ###########################################################################

    import numpy as np

    # 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.

            f.extend([0.,0.,0.])

        else:
            fi = f[-1]

            fip1 = min( 3.*g[i] , 3.*g[i+1] , np.sqrt(g[i+1]*g[i]) )

                # 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)

        # add next interval of interpolated (sub-)grid values
        f.append(fi1)
        f.append(fi2)
        f.append(fip1)

    return f