%********************************************************************
% This is a pure calculation (no plotting) script: EOF of sea level 
% anomalies along the equator for the period of 1961-2002 from the
% ERA40 exp. and NCEP exp. in Greatbatch et al (2018)
%
% Input:  eta3d_mm_fit.mat - 1961-2002 model-computed sea level
%	  Nino3_de.mat - 1961-2002 detrended Nino3 index 
%	  EMI_de.mat - 1961-2002 detrended EMI index
% Output: [input sea level file name]_EOF_eq.mat - EOF
%
% Written by Xiaoting Zhu, 2017 (xzhu@geomar.de)
%******************************************************************

%%--- BEGIN: global initialization ----------------------
clear all
close all
clc

crfolder = pwd  % current folder: script folder
rtfolder = [crfolder '/../'];  % root folder
dtfolder = [rtfolder 'Greatbatch2018_data/'];
subfolders = {'ERA40_exp/' 'NCEP_exp/'};

addpath([rtfolder 'mtools/EOF_analysis'])

%%--- Data file list ------------------------------------------
fname = {[dtfolder subfolders{1} 'eta3d_mm_fit.mat'],...
         [dtfolder subfolders{2} 'eta3d_mm_fit.mat']};

%%--- Loop: Read data, conduct EOF and save -----------------------------
for nf = 1:length(fname)
 load(fname{nf});
 time = time/24/3600+datenum('1957-09-01');
 datestr(time)

 %--- detrend and substract time mean
 etade = nan(size(eta));
 for i = 1:341
   for j = 1:61
    reta = squeeze(eta(i,j,:));
    retade = detrend(reta,'linear');
    etade(i,j,:) = retade - mean(retade);
   end
 end
 eta_de = squeeze(etade(:,24,:) + etade(:,25,:))/2.; % equator

 [n1,n3] = size(eta_de);
 ocean_msk = ~isnan(eta_de);
 X = eta_de(ocean_msk);

 n1 = size(X,1)/n3;

 X = reshape(X,[n1 n3]);

 %--- EOF: Xm*n m locations, n time, reverse of time-space
 [eigvet,L,pc] = EOF_analysis(X);

 eigval=100*L/sum(L);

 % modify sign and magnitude
 load([dtfolder subfolders{1} 'Nino3_de.mat']);
 [r,p]=corrcoef(squeeze(pc(1,:)),Nino3_de);
 sr = sign(r(1,2));
 pci=squeeze(pc(1,:));
 std(pci)
 eigvet(:,1) = sr*eigvet(:,1)*std(pci);
 pc(1,:) = sr*pci/std(pci);

 load([dtfolder subfolders{1} 'EMI_de.mat']);
 [r,p]=corrcoef(squeeze(pc(2,:)),EMI_de);
 sr = sign(r(1,2));
 pci=squeeze(pc(2,:));
 std(pci)
 eigvet(:,2) = sr*eigvet(:,2)*std(pci);
 pc(2,:) = sr*pci/std(pci);
 
%***********test*************
pcn=1;
eof=eigvet(:,pcn)*pc(pcn,:);
eof=eigvet*pc+2;
test=nan(size(eta_de));
test(ocean_msk)=eof;
recneof=test(301,:);
org=eta_de(301,:);
if nf==2
figure(1)
subplot(211)
plot([1:size(eof,2)],pc(pcn,:)/std(pc(pcn,:)),'k',[1:size(eof,2)],Nino3_de,'r',[1:size(eof,2)],EMI_de,'b')
title(['red Nino3, black PC' num2str(pcn) ', blue EMI_de'])
subplot(212)
plot([1:size(eof,2)],recneof,'r',[1:size(eof,2)],org,'k')
title(['red EOF' num2str(pcn) ', black total'])
title('red all EOF, black total')
cd(rtfolder)
saveas(gca,'test','epsc2')
end
%********************************
 
 % ocean mask
 eta2 = squeeze(eta_de(:,1));
 ocean_msk2 = ~isnan(eta2);
 num=2;

 %--- Save data 
 aa = fname{nf};
 save([aa(1:end-4) '_EOF_eq.mat'],'eigvet','pc','eigval','time','num','ocean_msk2') 

end
 
