What happened: With the latest dask release 2021.4.0, there seems to be a regression in memory management that has broken some of my standard climate science workflows like the simple anomaly calculation below. Rather than intelligently handling chunks for independent time slices, the code below will now load the entire ~30GB x array into memory before writing to zarr.

What you expected to happen: Dask would intelligently manage chunks and not fill up memory. This works fine in 2021.3.0.

Minimal Complete Verifiable Example: Generate a synthetic dataset with time/lat/lon variable and associated climatology stored to disk, then calculate the anomaly: ```python import xarray as xr import pandas as pd import numpy as np import dask.array as da

dates = pd.date_range('1980-01-01', '2019-12-31', freq='D') ds = xr.Dataset( data_vars = { 'x':( ('time', 'lat', 'lon'), da.random.random(size=(dates.size, 360, 720), chunks=(1, -1, -1))), 'clim':( ('dayofyear', 'lat', 'lon'), da.random.random(size=(366, 360, 720), chunks=(1, -1, -1))), }, coords = { 'time': dates, 'dayofyear': np.arange(1, 367, 1), 'lat': np.arange(-90, 90, .5), 'lon': np.arange(-180, 180, .5), } )

My original use case was pulling this data from disk, but it doesn't actually seem to matter

ds.to_zarr('test-data', mode='w') ds = xr.open_zarr('test-data') ds['anom'] = ds.x.groupby('time.dayofyear') - ds.clim ds[['anom']].to_zarr('test-anom', mode='w') ```

Anything else we need to know?: Distributed vs local scheduler and file backend e.g. zarr vs netcdf don't seem to affect this.

Dask graphs look the same for both 2021.3.0: and 2021.4.0:

Environment: