% Extract monthly averages
cd ~dmenemen/llc_4320/regions/Eq140W/matlab
nx=1; ny=832; nz=39; nt=13;
suf=['C_' int2str(nx) 'x' int2str(ny)];
lon=readbin(['../grid/X' suf],[nx ny]); lon=lon(1);
lat=readbin(['../grid/Y' suf],[nx ny]);
load ../grid/thk90
dpt=dpt90(1:nz);
ts=10368:144:1492992;
dte=datenum(2011,9,10)+ts*25/60/60/24;

for fld={'Salt','Theta','U','V','W'}
    pnm=['../' fld{1} '/'];
    if fld{1}=='V';
        suf=['_' fld{1} '_12385.7504.1_' int2str(nx) ...
             '.' int2str(ny) '.' int2str(nz) '_Neg'];
    else
        suf=['_' fld{1} '_12385.7505.1_' int2str(nx) ...
             '.' int2str(ny) '.' int2str(nz)];
    end
    eval([fld{1} '=zeros(ny,nz,nt);'])
    tme=zeros(nt,1);
    for m=10:22, mydisp(m)
        dte1=datenum(2011,m,1);
        dte2=datenum(2011,m+1,1);
        tme(m-9)=(dte1+dte2)/2;
        it=find(dte>=dte1&dte<=dte2);
        tmp=zeros(ny,nz);
        for t=1:length(it)
            fnm=[pnm myint2str(ts(it(t)),10) suf];
            tmp=tmp+readbin(fnm,[ny nz]);
        end
        eval([fld{1} '(:,:,m-9)=tmp/length(it);'])
    end
    eval(['save ' fld{1} 'monthly ' fld{1} ' lon lat dpt tme'])
    eval(['clear ' fld{1}])
end