function[Am,l,EOF]=eofv_grads(undef,perc,neof,norm,vin)
%-------------------------------------------------------------------
%[Am,l,EOF1,EOF2,...]=eofc_sm(undef,dim_x_y_t,perc,U1,U2,...)
%dim_x_y_t =[x_dim y_dim t_dim]
%U1(x,y,t), U2() <-- input fields
%perc <-- use only those grid points where data coverage exceeds perc[%]
%EOF1(x,y,#)
%EOF2(x,y,#)
%Am(t,#)
%l(#) eigenvalues
%-------------------------------------------------------------------
% Copyright (C) 2001 David M. Kaplan
% Licence: GPL (Gnu Public License)
% http://erizo.ucdavis.edu/~dmk/notes/EOFs/EOFs.html
% modified by senya@atmos.umd.edu for  multivariate  EOF
% v.1.2 tested with Matlab7.2 (2006a) 
%-------------------------------------------------------------------

[nx,ny,nz,nt,nvars]=size(vin);
    
if perc>1
    perc=perc/100;
end

vin=permute(vin,[4 1 2 3 5]);
vin=reshape(vin,[nt,nx*ny*nz,nvars]);


%create sparse matrix
validpoints=ones(1,nx*ny*nz);
dU=[];
%num=zeros(nx*ny,nvars);
for ivar=1:nvars
if isnan(undef)
S=sparse([~isnan(vin(:,:,ivar))].*unmaskout(vin(:,:,ivar),-9999));
%S=sparse([~isnan(varargin{ivar})].*unmaskout(varargin{ivar},-9999));
else
S=sparse([single(vin(:,:,ivar))~=single(undef)].*vin(:,:,ivar));
%S=sparse([single(varargin{ivar})~=single(undef)].*varargin{ivar});
end
SS=spones(S);

Um=sum(S)./sum(SS);
validpoints=validpoints & [sum(SS)>perc*max(sum(SS))];
dU(:,:,ivar)=(S-ones(nt,1)*Um).*SS;
end
ind=find(validpoints);
dU=dU(:,ind,:);





if norm==1
 for i=1:nvars
   norms(i)=std2(dU(:,:,i));
disp(['Variance' num2str(i) '='  num2str(norms(i))]);
   dU(:,:,i)=dU(:,:,i)/norms(i);
 end
    
end%if

dU=reshape(dU,nt,nvars*length(ind));
%
%totvar=std2(dU)^2; the same as below, but needs toolbox 
totvar=sum(sum(dU.^2))/(prod(size(dU))-1);
disp(['Total variance=' num2str(totvar)]);
clear  validpoints S SS

%mean zonal current at each point


clear Um 
[C,lambda,E]=svds(dU,neof);
l=diag(lambda.*lambda)/(nt);
Am=C*lambda;
explvar=100/totvar*l/(size(E,1));


for i=1:length(l)
E(:,i)=E(:,i)*sqrt(l(i));
Am(:,i)=Am(:,i)/sqrt(l(i));
end

for i=1:nvars
eof=NaN*ones(nx*ny*nz,size(E,2));
eof(ind,:)=E(1+(i-1)*length(ind):i*length(ind),:);
eof=reshape(eof,[nx,ny,nz,size(E,2)]);
if norm==1
  eof= eof*norms(i);
end %if
EOF(:,:,:,:,i)=eof;
end

Am(abs(Am)<1e-15)=NaN;
l=[l explvar];
