function [V_eof,PC,per,lamd]=EOF(A,kmod)
% [V_eof,PC,per,lamd]=EOF(A,kmod)
% A should be monthly field data (not anomaly) in [grid1,grid2,time] format
% kmode is the number of modes you'd like to have
Q_net=A;
clear A
[nx,ny,nt]=size(Q_net);
np=nx*ny;
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% reshape the data like [nx*ny,nt]
var=reshape(Q_net,np,nt);
% nn=mod(nt,12);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%calculate the anomaly of var
% for m=1:12;
% mon_mean(:,m)=mean(var(:,m:12:nt),2);
% if(m<= nn)
% var(:,m:12:nt)=var(:,m:12:nt)-mon_mean(:,m)*ones(1,floor(nt/12)+1);
% else
% var(:,m:12:nt)=var(:,m:12:nt)-mon_mean(:,m)*ones(1,floor(nt/12));
% end
% end
% for i=1:np
% var(i,:)=var(i,:)-mean(var(i,:));
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EOF analysis : V_eof , PC, per
%%%%% remove the point has NaN value from the data
%%%%% land is an index which is 1 if there is NaN means land
var_eof=var;
land(1:np)=0;
for i=1:np
index=find(isnan(var(i,:)));
if (length(index)~=0); % 对某个空间点np，在nt的时间序列中，只要有一个时间点上有陆地值或不合理的数据值，则所有时间点上，此空间点都算作陆地值(即无值点)
land(i)=1; %%%%%%%% this means here is land
end
end
if(isempty(find(land==1))==1)
else
var_eof(find(land==1),:)=[]; %%%% delete the point where there is NaN
end
%for i=1:k-1, var_eof(i,:)=var_eof(i,:)-mean(var_eof(i,:)); end;
%%%%%%%%
%%%%% var_eof is the data to EOF
% kmod=4; %%%%%%%%%%%%%%%%%% the number of mode wanted
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[V1,E]=eig(var_eof'*var_eof); % 空间点大于时间点的时候为了方便计算，用时空转换的方法求特征值和特征向量
E=fliplr(flipud(E));
V2=fliplr(var_eof*V1);
for i=1:kmod;
V(:,i)=V2(:,i) / sqrt(E(i,i));
PC(:,i)=(V(:,i)'* var_eof)';
lamd=diag(E); % 取主对角线上的元素
per(i)=E(i,i) / sum(lamd);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% put V back as land to Vl (Vland) ,V_eof(nx,ny,kmod)
Vl(1:np,1:kmod)=NaN; k=1;
for i=1:np;
if land(i)==0 ; Vl(i,:)=V(k,:); k=k+1; end;
end;
V_eof=reshape(Vl,nx,ny,kmod);
%
