def fa(arr, t, dist='norm', mode='high'):
    """Return the value corresponding to the given return period.

    Parameters
    ----------
    arr : xarray.DataArray
      Maximized/minimized input data with a `time` dimension.
    t : int or sequence
      Return period. The period depends on the resolution of the input data. If the input array's resolution is
      yearly, then the return period is in years.
    dist : str
      Name of the univariate distribution, such as beta, expon, genextreme, gamma, gumbel_r, lognorm, norm
      (see scipy.stats).
    mode : {'min', 'max}
      Whether we are looking for a probability of exceedance (max) or a probability of non-exceedance (min).

    Returns
    -------
    xarray.DataArray
      An array of values with a 1/t probability of exceedance (if mode=='max').
    """
    t = np.atleast_1d(t)

    # Get the distribution
    dc = get_dist(dist)

    # Fit the parameters of the distribution
    p = fit(arr, dist)

    # Create a lambda function to facilitate passing arguments to dask. There is probably a better way to do this.
    if mode in ['max', 'high']:
        def func(x):
            return dc.isf(1./t, *x)
    elif mode in ['min', 'low']:
        def func(x):
            return dc.ppf(1./t, *x)
    else:
        raise ValueError("mode `{}` should be either 'max' or 'min'".format(mode))

    data = dask.array.apply_along_axis(func, p.get_axis_num('dparams'), p)

    # Create coordinate for the return periods
    coords = dict(p.coords.items())
    coords.pop('dparams')
    coords['return_period'] = t

    # Create dimensions
    dims = list(p.dims)
    dims.remove('dparams')
    dims.insert(0, u'return_period')

    # TODO: add time and time_bnds coordinates (Low will work on this)
    # time.attrs['climatology'] = 'climatology_bounds'
    # coords['time'] =
    # coords['climatology_bounds'] =

    out = xr.DataArray(data=data, coords=coords, dims=dims)
    out.attrs = p.attrs
    out.attrs['standard_name'] = '{0} quantiles'.format(dist)
    out.attrs['long_name'] = '{0} return period values for {1}'.format(dist, getattr(arr, 'standard_name', ''))
    out.attrs['cell_methods'] = (out.attrs.get('cell_methods', '') + ' dparams: ppf').strip()
    out.attrs['units'] = arr.attrs.get('units', '')
    out.attrs['mode'] = mode
    out.attrs['history'] = out.attrs.get('history', '') + "Compute values corresponding to return periods."

    return out
Return the value corresponding to the given return period.

    Parameters
    ----------
    arr : xarray.DataArray
      Maximized/minimized input data with a `time` dimension.
    t : int or sequence
      Return period. The period depends on the resolution of the input data. If the input array's resolution is
      yearly, then the return period is in years.
    dist : str
      Name of the univariate distribution, such as beta, expon, genextreme, gamma, gumbel_r, lognorm, norm
      (see scipy.stats).
    mode : {'min', 'max}
      Whether we are looking for a probability of exceedance (max) or a probability of non-exceedance (min).

    Returns
    -------
    xarray.DataArray
      An array of values with a 1/t probability of exceedance (if mode=='max').