Create an xarray datset of ensemble of climate simulation from a list of netcdf files. Input data is
    concatenated along a newly created data dimension ('realization')

    Returns a xarray dataset object containing input data from the list of netcdf files concatenated along
    a new dimension (name:'realization'). In the case where input files have unequal time dimensions output
    ensemble dataset is created for overlapping time-steps common to all input files

    Parameters
    ----------
    ncfiles : sequence
      List of netcdf file paths. If mf_flag is true ncfiles should be a list of lists where
    each sublist contains input .nc files of a multifile dataset

    mf_flag : Boolean . If true climate simulations are treated as multifile datasets before concatenation

    Returns
    -------
    xarray dataset containing concatenated data from all input files

    Notes
    -----
    Input netcdf files require equal spatial dimension size (e.g. lon, lat dimensions)
    If input data contains multiple cftime calendar types they must be at monthly or coarser frequency

    Examples
    --------
    >>> from xclim import utils
    >>> import glob
    >>> ncfiles = glob.glob('/*.nc')
    >>> ens = utils.create_ensemble(ncfiles)
    >>> print(ens)
    Using multifile datasets:
    simulation 1 is a list of .nc files (e.g. separated by time)
    >>> ncfiles = glob.glob('dir/*.nc')
    simulation 2 is also a list of .nc files
    >>> ens = utils.create_ensemble(ncfiles)

def create_ensemble(ncfiles, mf_flag=False):
    """Create an xarray datset of ensemble of climate simulation from a list of netcdf files. Input data is
    concatenated along a newly created data dimension ('realization')

    Returns a xarray dataset object containing input data from the list of netcdf files concatenated along
    a new dimension (name:'realization'). In the case where input files have unequal time dimensions output
    ensemble dataset is created for overlapping time-steps common to all input files

    Parameters
    ----------
    ncfiles : sequence
      List of netcdf file paths. If mf_flag is true ncfiles should be a list of lists where
    each sublist contains input .nc files of a multifile dataset

    mf_flag : Boolean . If true climate simulations are treated as multifile datasets before concatenation

    Returns
    -------
    xarray dataset containing concatenated data from all input files

    Notes
    -----
    Input netcdf files require equal spatial dimension size (e.g. lon, lat dimensions)
    If input data contains multiple cftime calendar types they must be at monthly or coarser frequency

    Examples
    --------
    >>> from xclim import utils
    >>> import glob
    >>> ncfiles = glob.glob('/*.nc')
    >>> ens = utils.create_ensemble(ncfiles)
    >>> print(ens)
    Using multifile datasets:
    simulation 1 is a list of .nc files (e.g. separated by time)
    >>> ncfiles = glob.glob('dir/*.nc')
    simulation 2 is also a list of .nc files
    >>> ens = utils.create_ensemble(ncfiles)
    """
    dim = 'realization'
    ds1 = []
    start_end_flag = True
    print('finding common time-steps')
    for n in ncfiles:
        if mf_flag:
            ds = xr.open_mfdataset(n, concat_dim='time', decode_times=False, chunks={'time': 10})
            ds['time'] = xr.open_mfdataset(n).time
        else:
            ds = xr.open_dataset(n, decode_times=False)
            ds['time'] = xr.decode_cf(ds).time
        # get times - use common
        time1 = pd.to_datetime({'year': ds.time.dt.year, 'month': ds.time.dt.month, 'day': ds.time.dt.day})
        if start_end_flag:
            start1 = time1.values[0]
            end1 = time1.values[-1]
            start_end_flag = False
        if time1.values.min() > start1:
            start1 = time1.values.min()
        if time1.values.max() < end1:
            end1 = time1.values.max()

    for n in ncfiles:
        print('accessing file ', ncfiles.index(n) + 1, ' of ', len(ncfiles))
        if mf_flag:
            ds = xr.open_mfdataset(n, concat_dim='time', decode_times=False, chunks={'time': 10})
            ds['time'] = xr.open_mfdataset(n).time
        else:
            ds = xr.open_dataset(n, decode_times=False, chunks={'time': 10})
            ds['time'] = xr.decode_cf(ds).time

        ds['time'].values = pd.to_datetime({'year': ds.time.dt.year, 'month': ds.time.dt.month, 'day': ds.time.dt.day})

        ds = ds.where((ds.time >= start1) & (ds.time <= end1), drop=True)

        ds1.append(ds.drop('time'))
    print('concatenating files : adding dimension ', dim, )
    ens = xr.concat(ds1, dim=dim)
    # assign time coords
    ens = ens.assign_coords(time=ds.time.values)
    return ens