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);

vrhom=zeros(141,361,55,12);


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

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

filename

eval(['load ' filename])

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

vrhom(:,:,:,imon)=vrhom(:,:,:,imon)+nanmean(vv,4);

end
end

vrhom=vrhom/length(year);

filename=[datadir 'vrhom'];

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