def interpolate_variable(
    var: Variable,
    indexes_coords: SourceDest,
    *,
    method: InterpOptions,
    kwargs: dict[str, Any],
) -> Variable:
    """core routine that returns the interpolated variable."""
    if not indexes_coords:
        return var.copy()

    if len(indexes_coords) == 1:
        func, kwargs = _get_interpolator(method, vectorizeable_only=True, **kwargs)
    else:
        func, kwargs = _get_interpolator_nd(method, **kwargs)

    in_coords, result_coords = zip(*(v for v in indexes_coords.values()), strict=True)

    # input coordinates along which we are interpolation are core dimensions
    # the corresponding output coordinates may or may not have the same name,
    # so `all_in_core_dims` is also `exclude_dims`
    all_in_core_dims = set(indexes_coords)

    result_dims = OrderedSet(itertools.chain(*(_.dims for _ in result_coords)))
    result_sizes = ChainMap(*(_.sizes for _ in result_coords))

    # any dimensions on the output that are present on the input, but are not being
    # interpolated along are dimensions along which we automatically vectorize.
    # Consider the problem in https://github.com/pydata/xarray/issues/6799#issuecomment-2474126217
    # In the following, dimension names are listed out in [].
    # # da[time, q, lat, lon].interp(q=bar[lat,lon]). Here `lat`, `lon`
    # are input dimensions, present on the output, but are not the coordinates
    # we are explicitly interpolating. These are the dimensions along which we vectorize.
    # `q` is the only input core dimensions, and changes size (disappears)
    # so it is in exclude_dims.
    vectorize_dims = (result_dims - all_in_core_dims) & set(var.dims)

    # remove any output broadcast dimensions from the list of core dimensions
    output_core_dims = tuple(d for d in result_dims if d not in vectorize_dims)
    input_core_dims = (
        # all coordinates on the input that we interpolate along
        [tuple(indexes_coords)]
        # the input coordinates are always 1D at the moment, so we just need to list out their names
        + [tuple(_.dims) for _ in in_coords]
        # The last set of inputs are the coordinates we are interpolating to.
        + [
            tuple(d for d in coord.dims if d not in vectorize_dims)
            for coord in result_coords
        ]
    )
    output_sizes = {k: result_sizes[k] for k in output_core_dims}

    # scipy.interpolate.interp1d always forces to float.
    dtype = float if not issubclass(var.dtype.type, np.inexact) else var.dtype
    result = apply_ufunc(
        _interpnd,
        var,
        *in_coords,
        *result_coords,
        input_core_dims=input_core_dims,
        output_core_dims=[output_core_dims],
        exclude_dims=all_in_core_dims,
        dask="parallelized",
        kwargs=dict(
            interp_func=func,
            interp_kwargs=kwargs,
            # we leave broadcasting up to dask if possible
            # but we need broadcasted values in _interpnd, so propagate that
            # context (dimension names), and broadcast there
            # This would be unnecessary if we could tell apply_ufunc
            # to insert size-1 broadcast dimensions
            result_coord_core_dims=input_core_dims[-len(result_coords) :],
        ),
        # TODO: deprecate and have the user rechunk themselves
        dask_gufunc_kwargs=dict(output_sizes=output_sizes, allow_rechunk=True),
        output_dtypes=[dtype],
        vectorize=bool(vectorize_dims),
        keep_attrs=True,
    )
    return result
core routine that returns the interpolated variable.
xarray/core/missing.py