def interpolate_1d(x, xp, *args, **kwargs):
    r
    fill_value = kwargs.pop('fill_value', np.nan)
    axis = kwargs.pop('axis', 0)
    x = np.asanyarray(x).reshape(-1)
    ndim = xp.ndim
    sort_args = np.argsort(xp, axis=axis)
    sort_x = np.argsort(x)
    sorter = broadcast_indices(xp, sort_args, ndim, axis)
    xp = xp[sorter]
    variables = [arr[sorter] for arr in args]
    x_array = x[sort_x]
    expand = [np.newaxis] * ndim
    expand[axis] = slice(None)
    x_array = x_array[tuple(expand)]
    minv = np.apply_along_axis(np.searchsorted, axis, xp, x[sort_x])
    minv2 = np.copy(minv)
    if ((np.max(minv) == xp.shape[axis]) or (np.min(minv) == 0)) and fill_value is None:
        raise ValueError('Interpolation point out of data bounds encountered')
    if np.max(minv) == xp.shape[axis]:
        warnings.warn('Interpolation point out of data bounds encountered')
        minv2[minv == xp.shape[axis]] = xp.shape[axis] - 1
    if np.min(minv) == 0:
        minv2[minv == 0] = 1
    above = broadcast_indices(xp, minv2, ndim, axis)
    below = broadcast_indices(xp, minv2 - 1, ndim, axis)
    if np.any(x_array < xp[below]):
        warnings.warn('Interpolation point out of data bounds encountered')
    ret = []
    for var in variables:
        var_interp = var[below] + (var[above] - var[below]) * ((x_array - xp[below])
                                                               / (xp[above] - xp[below]))
        var_interp[minv == xp.shape[axis]] = fill_value
        var_interp[x_array < xp[below]] = fill_value
        if x[0] > x[-1]:
            var_interp = np.swapaxes(np.swapaxes(var_interp, 0, axis)[::-1], 0, axis)
        ret.append(var_interp)
    if len(ret) == 1:
        return ret[0]
    else:
        return ret