% Extract monthly averages
cd ~dmenemen/llc_2160/regions/Eq140Wbox/matlab
nx=120; ny=364; nz=36; nt=24;
suf=['C_' int2str(nx) 'x' int2str(ny)];
lon=readbin(['../grid/X' suf],[nx ny]); lon=lon(:,1);
lat=readbin(['../grid/Y' suf],[nx ny]); lat=lat(1,:);
load ~dmenemen/llc_4320/grid/thk90.mat
dpt=dpt90(1:nz);
ts=92160:80:1586400;
dte=datenum(2011,1,17)+ts*45/60/60/24;

for fld={'Salt','Theta','U','V','W'}
    pnm=['../' fld{1} '/'];
    if fld{1}=='V';
        suf=['_' fld{1} '_6133.3804.1_' int2str(nx) ...
             '.' int2str(ny) '.' int2str(nz)];
    else
        suf=['_' fld{1} '_6133.3805.1_' int2str(nx) ...
             '.' int2str(ny) '.' int2str(nz)];
    end
    eval([fld{1} '=zeros(nx,ny,nz,nt);'])
    tme=zeros(nt,1);
    for m=4:27, disp(m)
        dte1=datenum(2011,m,1);
        dte2=datenum(2011,m+1,1);
        tme(m-3)=(dte1+dte2)/2;
        it=find(dte>=dte1&dte<=dte2);
        tmp=zeros(nx,ny,nz);
        for t=1:length(it)
            fnm=[pnm myint2str(ts(it(t)),10) suf];
            tmp=tmp+readbin(fnm,[nx ny nz]);
        end
        eval([fld{1} '(:,:,:,m-3)=tmp/length(it);'])
    end
    eval(['save ' fld{1} '_monthly ' fld{1} ' lon lat dpt tme'])
    eval(['clear ' fld{1}])
end