concat_dim getting added to *all* variables of multifile datasets

Code Sample

Using the following example data set:
example_jan.nc

#!/usr/bin/env python3

import xarray
ds = xarray.open_mfdataset('example_jan.nc', concat_dim='Time')
print(ds)

The result from xarray 0.10.2 (and all previous various xarray versions we've worked with):

Dimensions:                                                      (Time: 1, nOceanRegions: 7, nOceanRegionsTmp: 7, nVertLevels: 100)
Dimensions without coordinates: Time, nOceanRegions, nOceanRegionsTmp, nVertLevels
Data variables:
    time_avg_avgValueWithinOceanLayerRegion_avgLayerTemperature  (Time, nOceanRegionsTmp, nVertLevels) float64 dask.array<shape=(1, 7, 100), chunksize=(1, 7, 100)>
    time_avg_avgValueWithinOceanRegion_avgSurfaceTemperature     (Time, nOceanRegions) float64 dask.array<shape=(1, 7), chunksize=(1, 7)>
    time_avg_daysSinceStartOfSim                                 (Time) timedelta64[ns] dask.array<shape=(1,), chunksize=(1,)>
    xtime_end                                                    (Time) |S64 dask.array<shape=(1,), chunksize=(1,)>
    xtime_start                                                  (Time) |S64 dask.array<shape=(1,), chunksize=(1,)>
    refBottomDepth                                               (nVertLevels) float64 dask.array<shape=(100,), chunksize=(100,)>
Attributes:
    history:  Tue Dec  6 04:49:14 2016: ncatted -O -a ,global,d,, acme_alaph7...
    NCO:      "4.6.2"

The results with xarray 0.10.3:

<xarray.Dataset>
Dimensions:                                                      (Time: 1, nOceanRegions: 7, nOceanRegionsTmp: 7, nVertLevels: 100)
Dimensions without coordinates: Time, nOceanRegions, nOceanRegionsTmp, nVertLevels
Data variables:
    time_avg_avgValueWithinOceanLayerRegion_avgLayerTemperature  (Time, nOceanRegionsTmp, nVertLevels) float64 dask.array<shape=(1, 7, 100), chunksize=(1, 7, 100)>
    time_avg_avgValueWithinOceanRegion_avgSurfaceTemperature     (Time, nOceanRegions) float64 dask.array<shape=(1, 7), chunksize=(1, 7)>
    time_avg_daysSinceStartOfSim                                 (Time) timedelta64[ns] dask.array<shape=(1,), chunksize=(1,)>
    xtime_end                                                    (Time) |S64 dask.array<shape=(1,), chunksize=(1,)>
    xtime_start                                                  (Time) |S64 dask.array<shape=(1,), chunksize=(1,)>
    refBottomDepth                                               (Time, nVertLevels) float64 dask.array<shape=(1, 100), chunksize=(1, 100)>
Attributes:
    history:  Tue Dec  6 04:49:14 2016: ncatted -O -a ,global,d,, acme_alaph7...
    NCO:      "4.6.2"

Problem description

The expected behavior for us was that refBottomDepth should not have Time as a dimension. It does not vary with time and does not have a Time dimension in the input data set.

It seems like https://github.com/pydata/xarray/issues/1988 and https://github.com/pydata/xarray/pull/2048 were intended to address cases where the concat_dim was not yet present in the input files. But in cases where concat_dim is already in the input files, it seems like only those fields that include this dimensions should be concatenated and other fields should remain free of concat_dim. Part of the problem for us is that the number of dimensions of some of our variables change depending on which xarray version is being used.

Expected Output

That for 0.10.2 (see above)

Output of xr.show_versions()

Reproduce the behavior with xarray.open_mfdataset('example_jan.nc', concat_dim='Time') using the linked example_jan.nc dataset, and compare the refBottomDepth dimensions across xarray 0.10.2 and 0.10.3. Trace the multifile concatenation path to determine why concat_dim is added to variables that lack it; done means refBottomDepth remains (nVertLevels) while variables containing Time still concatenate correctly.