def _dataset_concat(
    datasets: Iterable[T_Dataset],
    dim: str | T_Variable | T_DataArray | pd.Index,
    data_vars: T_DataVars,
    coords: str | list[str],
    compat: CompatOptions,
    positions: Iterable[Iterable[int]] | None,
    fill_value: Any = dtypes.NA,
    join: JoinOptions = "outer",
    combine_attrs: CombineAttrsOptions = "override",
    create_index_for_new_dim: bool = True,
) -> T_Dataset:
    """
    Concatenate a sequence of datasets along a new or existing dimension
    """
    from xarray.core.dataarray import DataArray
    from xarray.core.dataset import Dataset

    datasets = list(datasets)

    if not all(isinstance(dataset, Dataset) for dataset in datasets):
        raise TypeError(
            "The elements in the input list need to be either all 'Dataset's or all 'DataArray's"
        )

    if isinstance(dim, DataArray):
        dim_var = dim.variable
    elif isinstance(dim, Variable):
        dim_var = dim
    else:
        dim_var = None

    dim_name, index = _calc_concat_dim_index(dim)

    # Make sure we're working on a copy (we'll be loading variables)
    datasets = [ds.copy() for ds in datasets]
    datasets = list(
        align(
            *datasets, join=join, copy=False, exclude=[dim_name], fill_value=fill_value
        )
    )

    dim_coords, dims_sizes, coord_names, data_names, vars_order = _parse_datasets(
        datasets
    )
    dim_names = set(dim_coords)

    both_data_and_coords = coord_names & data_names
    if both_data_and_coords:
        raise ValueError(
            f"{both_data_and_coords!r} is a coordinate in some datasets but not others."
        )
    # we don't want the concat dimension in the result dataset yet
    dim_coords.pop(dim_name, None)
    dims_sizes.pop(dim_name, None)

    # case where concat dimension is a coordinate or data_var but not a dimension
    if (
        dim_name in coord_names or dim_name in data_names
    ) and dim_name not in dim_names:
        datasets = [
            ds.expand_dims(dim_name, create_index_for_new_dim=create_index_for_new_dim)
            for ds in datasets
        ]

    # determine which variables to concatenate
    concat_over, equals, concat_dim_lengths, skip_merge = _calc_concat_over(
        datasets, dim_name, dim_names, data_vars, coords, compat
    )

    # determine which variables to merge, and then merge them according to compat
    variables_to_merge = (coord_names | data_names) - concat_over - skip_merge

    result_vars = {}
    result_indexes = {}

    if variables_to_merge:
        grouped = {
            k: v
            for k, v in collect_variables_and_indexes(datasets).items()
            if k in variables_to_merge
        }
        merged_vars, merged_indexes = merge_collected(
            grouped, compat=compat, equals=equals
        )
        result_vars.update(merged_vars)
        result_indexes.update(merged_indexes)

    result_vars.update(dim_coords)

    # assign attrs and encoding from first dataset
    result_attrs = merge_attrs([ds.attrs for ds in datasets], combine_attrs)
    result_encoding = datasets[0].encoding

    # check that global attributes are fixed across all datasets if necessary
    if compat == "identical":
        for ds in datasets[1:]:
            if not utils.dict_equiv(ds.attrs, result_attrs):
                raise ValueError("Dataset global attributes not equal.")

    # we've already verified everything is consistent; now, calculate
    # shared dimension sizes so we can expand the necessary variables
    def ensure_common_dims(vars, concat_dim_lengths):
        # ensure each variable with the given name shares the same
        # dimensions and the same shape for all of them except along the
        # concat dimension
        common_dims = tuple(utils.OrderedSet(d for v in vars for d in v.dims))
        if dim_name not in common_dims:
            common_dims = (dim_name,) + common_dims
        for var, dim_len in zip(vars, concat_dim_lengths, strict=True):
            if var.dims != common_dims:
                common_shape = tuple(dims_sizes.get(d, dim_len) for d in common_dims)
                var = var.set_dims(common_dims, common_shape)
            yield var

    # get the indexes to concatenate together, create a PandasIndex
    # for any scalar coordinate variable found with ``name`` matching ``dim``.
    # TODO: depreciate concat a mix of scalar and dimensional indexed coordinates?
    # TODO: (benbovy - explicit indexes): check index types and/or coordinates
    # of all datasets?
    def get_indexes(name):
        for ds in datasets:
            if name in ds._indexes:
                yield ds._indexes[name]
            elif name == dim_name:
                var = ds._variables[name]
                if not var.dims:
                    data = var.set_dims(dim_name).values
                    if create_index_for_new_dim:
                        yield PandasIndex(data, dim_name, coord_dtype=var.dtype)

    # create concatenation index, needed for later reindexing
    file_start_indexes = np.append(0, np.cumsum(concat_dim_lengths))
    concat_index = np.arange(file_start_indexes[-1])
    concat_index_size = concat_index.size
    variable_index_mask = np.ones(concat_index_size, dtype=bool)

    # stack up each variable and/or index to fill-out the dataset (in order)
    # n.b. this loop preserves variable order, needed for groupby.
    ndatasets = len(datasets)
    for name in vars_order:
        if name in concat_over and name not in result_indexes:
            variables = []
            # Initialize the mask to all True then set False if any name is missing in
            # the datasets:
            variable_index_mask.fill(True)
            var_concat_dim_length = []
            for i, ds in enumerate(datasets):
                if name in ds.variables:
                    variables.append(ds[name].variable)
                    var_concat_dim_length.append(concat_dim_lengths[i])
                else:
                    # raise if coordinate not in all datasets
                    if name in coord_names:
                        raise ValueError(
                            f"coordinate {name!r} not present in all datasets."
                        )

                    # Mask out the indexes without the name:
                    start = file_start_indexes[i]
                    end = file_start_indexes[i + 1]
                    variable_index_mask[slice(start, end)] = False

            variable_index = concat_index[variable_index_mask]
            vars = ensure_common_dims(variables, var_concat_dim_length)

            # Try to concatenate the indexes, concatenate the variables when no index
            # is found on all datasets.
            indexes: list[Index] = list(get_indexes(name))
            if indexes:
                if len(indexes) < ndatasets:
                    raise ValueError(
                        f"{name!r} must have either an index or no index in all datasets, "
                        f"found {len(indexes)}/{len(datasets)} datasets with an index."
                    )
                combined_idx = indexes[0].concat(indexes, dim_name, positions)
                if name in datasets[0]._indexes:
                    idx_vars = datasets[0].xindexes.get_all_coords(name)
                else:
                    # index created from a scalar coordinate
                    idx_vars = {name: datasets[0][name].variable}
                result_indexes.update(dict.fromkeys(idx_vars, combined_idx))
                combined_idx_vars = combined_idx.create_variables(idx_vars)
                for k, v in combined_idx_vars.items():
                    v.attrs = merge_attrs(
                        [ds.variables[k].attrs for ds in datasets],
                        combine_attrs=combine_attrs,
                    )
                    result_vars[k] = v
            else:
                combined_var = concat_vars(
                    vars, dim_name, positions, combine_attrs=combine_attrs
                )
                # reindex if variable is not present in all datasets
                if len(variable_index) < concat_index_size:
                    combined_var = reindex_variables(
                        variables={name: combined_var},
                        dim_pos_indexers={
                            dim_name: pd.Index(variable_index).get_indexer(concat_index)
                        },
                        fill_value=fill_value,
                    )[name]
                result_vars[name] = combined_var

        elif name in result_vars:
            # preserves original variable order
            result_vars[name] = result_vars.pop(name)

    absent_coord_names = coord_names - set(result_vars)
    if absent_coord_names:
        raise ValueError(
            f"Variables {absent_coord_names!r} are coordinates in some datasets but not others."
        )

    result_data_vars = {}
    coord_vars = {}
    for name, result_var in result_vars.items():
        if name in coord_names:
            coord_vars[name] = result_var
        else:
            result_data_vars[name] = result_var

    if index is not None:
        if dim_var is not None:
            index_vars = index.create_variables({dim_name: dim_var})
        else:
            index_vars = index.create_variables()

        coord_vars[dim_name] = index_vars[dim_name]
        result_indexes[dim_name] = index

    coords_obj = Coordinates(coord_vars, indexes=result_indexes)

    result = type(datasets[0])(result_data_vars, coords=coords_obj, attrs=result_attrs)
    result.encoding = result_encoding

    return result

    Concatenate a sequence of datasets along a new or existing dimension
    
xarray/structure/concat.py