'Time' dimensions not consistently saved when using 'to_netcdf'

What happened:
When using a for loop to subset 3 datasets, the second dataset only preserves 1 Time dimension while the first and third preserve all 4 Time dimensions. However, when subsetting the datasets individually (i.e. not using a for loop), the Time dimensions are all preserved for each file. I'm not sure if this is a bug or user error, but could possibly be related to #5549.

What you expected to happen:
Expected all three files to subset in the same capacity with same dimensions for variables.

Minimal Complete Verifiable Example:
Python code for the for loop to subset:

# Put your MCVE code here

wrf_path = '/Users/misi1684/wrf_python_tutorial/wrf_tutorial_data'
os.chdir(wrf_path)
wrf_file0 = glob.glob('wrf*')

toinclude = ['HGT', 'QVAPOR', 'QRAIN', 'QCLOUD', 'P', 'PB', 'PHB', 'PH', 'T', 'U', 'V', 'W']

for i in range(0,3):
	wrf_file = wrf_file0[i]
	wrf_out = xr.open_dataset(wrf_file, decode_times=False)
	# print(wrf_out)
	wrf_out = wrf_out[toinclude].to_netcdf('wrf_output_subset' + str(i) + '.nc')

Anything else we need to know?:
Incase you are not familiar with a wrfout file, it is the same as a netcdf file. Generally, I used netCDF Dataset to read in these types of files, but netCDF does not have the same subsetting capabilities as xarray. If you need more information or printouts from these datasets, please let me know!

Environment:

Start by running the provided Python loop with the three WRF/netCDF datasets and compare the dimensions in each generated file with individually subset files. Investigate the to_netcdf path involved in preserving dimensions. Done means all three subset files preserve the same expected Time dimensions.