%climatological eddy kinetic energy

exp='noed'

idir=['/Volumes/scratch/baedw/data/processed/' exp];

ufils=dir([idir '/usurf*']);
vfils=dir([idir '/vsurf*']);

  ncload([idir '/' ufils(1).name]);
  UU=zeros(length(yu_ocean),length(xu_ocean), 12);
  VV=UU;
  Uab=UU;
  PE=UU;
  EE=UU;

for k=[1:12]
 U=zeros(length(yu_ocean), length(xu_ocean), 1);
 V=zeros(length(yu_ocean), length(xu_ocean), 1);
 n=0;
 for i=1:length(ufils)
  display(['processing month ' num2str(k) '; loading ' ufils(i).name])
  ncload([idir '/' ufils(i).name]); 
  ncload([idir '/' vfils(i).name]); 
  u=squee(u); v=squee(v);
  e=0.5.*(u.^2+v.^2).^0.5; 
  I=round(linspace(1,length(time),13));
  for l=I(k):I(k+1)
   n=n+1;
   U(:,:,n)=squee(u(l,:,:));
   V(:,:,n)=squee(v(l,:,:));
  end
 end
   UU(:,:,k)=mean(U,3);
   VV(:,:,k)=mean(V,3);
   Uab(:,:,k)=mean(0.5.*(U.^2+V.^2).^0.5,3);
   PE(:,:,k)=(UU(:,:,k).^2+VV(:,:,k).^2).^0.5 ./ mean((U.^2+V.^2).^0.5,3);
   for i=1:length(xu_ocean)
    for j=1:length(yu_ocean);
     EE(j,i,k)=0.5.*mean( (squee(U(j,i,:))-squee(UU(j,i,k))).^2 + (squee(V(j,i,:))-squee(VV(j,i,k))).^2  );
    end
   end
end

I=find(UU<-1e10); UU(I)=NaN; EE(I)=NaN; VV(I)=NaN; PE(I)=NaN; Uab(I)=NaN;
x=xu_ocean;
y=yu_ocean;

save([idir '/surfmom.mat'],'x','y','UU','VV','EE','PE','Uab')
