clear

casenum=input('Case no: ');

mon = char('jan','feb','mar','apr','may','jun','jul','aug','sep','oct','nov','dec');
nday = [31 28 31 30 31 30 31 31 30 31 30 31];
year = [1999 2000 2001 2002];

datadir = sprintf('/export/t2d2/kelvin/iROAM/daily_output/case%d/',casenum);

hum=zeros(141,361,56,12);

for iyear=1:length(year)
for imon=1:12

filename=[datadir 'daily_hu_' mon(imon,:) '_' int2str(year(iyear))];

filename

eval(['load ' filename])

% change order to original (from DODS)
hu=permute(hu,[4,3,2,1]);

hum(:,:,:,imon)=hum(:,:,:,imon)+nanmean(hu,4);

end
end

hum=hum/length(year);

filename=[datadir 'hum'];

eval(['save ' filename ' hum'])
