It's hard to decipher your nested loops, but here's how I'd tackle it. I'll assume you have a big 3D matrix U whose dimensions correspond to (lon x lat x time) or (gridx x gridy x time) or something of the like. You probably also have a 3D matrix of the same size V for zonal motion, and I'll assume you have an array t and length(t) = size(U,3).
1. Get t into components of year, month, and day:
[year,month,day] = datevec(t);
2. The downsample_ts function computes monthly means, so we'll lie to downsample_ts and tell it that all December and February data happened in January.
month(month==12) = 1;
month(month==2) = 1;
Do the same for the other seasons, i.e.,
month(month==3) = 4;
month(month==5) = 5;
And when you do the same for the other seasons, you should get
unique(month)
= [1, 4, 7, 10];
Create a new fake t2 vector:
t2 = datenum([year,month,day]);
Then a 32 year record of seasonal means can be obtained by
[Umean,tmean] = downsample_ts(U,t,'monthly','mean');
Vmean = downsample_ts(V,t,'monthly','mean');
Dimensions of Umean and Vmean will be size(U,1) x size(U,2) x 128 because 4 seasons times 32 years is 128. All DJF data will be given by
Umean(:,:,tmean==1);
And if you want to get a linear trend map of DJF seasonal averages over 32 years you could do
U_djf_trend = trend(Umean(:,:,tmean==1),1,3);
