%******************************************************************
% This script is for Figure 9
%
% Input:  EOF_eq_mm_fit.mat - from ERA40 exp. and NCEP exp.
%	  EOF_eq_de.mat  - from AVISO and Zhu et al (2017)
%	  EP_[ERA...]_eastward_stress_modelgrid.nc - from ERA40 or ERA-Interim
%	  EP_ustr_anom_mon_mean_NCEP_modelgrid.nc - from NCEP
% Output: Figure9.eps
% Note:   After s0_nodetrend_fit_mm_to_aviso_Nmodes_readycoeff.m
%	        s2_detrend_EOF_eq_Zhu2017.m
%		s3_detrend_EOF_eq_Greatbatch2018.m
%
% Written by Xiaoting Zhu, 2017 (xzhu@geomar.de)
%******************************************************************

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

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

set(0,'defaultfigurecolor',[1 1 1])
addpath([rtfolder 'mtools/altmany-export_fig-9ac0917']);
addpath([rtfolder 'mtools/Zhu_mfunct'])

lont=[120:20:290];
lonstr={'120\circE' '140\circE' '160\circE' '180\circ'...
        '160\circW' '140\circW' '120\circW' '100\circW' '80\circW'};
    
%%--- Data file list --------------- 
fname = {[dtfolder subfolders{1} 'eta3d_mm_fit_EOF_eq.mat'],...
	 [dtfolder subfolders{2} 'eta3d_mm_fit_EOF_eq.mat'],...
	 [dtfolder subfolders{3} 'sla3d_EOF_eq.mat'],...
	 [dtfolder subfolders{4} 'eta3d_mm_fit_EOF_eq.mat']};

fnamew= {[dtfolder subfolders{1} 'EP_ERA40_anom_mon_eastward_stress_modelgrid.nc'],...
	 [dtfolder subfolders{2} 'EP_ustr_anom_mon_mean_NCEP_modelgrid.nc'],...
	 [dtfolder subfolders{4} 'EP_ERA_interim_anom_mon_eastward_stress_modelgrid.nc'],...
	 [dtfolder subfolders{4} 'EP_ERA_interim_anom_mon_eastward_stress_modelgrid.nc']}

%%--- Read data and plot ----------------
hl=nan(1,length(fname)); %1-4
linesty = {'r','r--','k','k--'};
source={'ERA-40 exp. 61-02',...
        'NCEP exp. 61-02',...
        'AVISO 93-14',...
        'ERA-I exp. 93-14'};

figure('position',[100 100 500 200])
lon = ncread([dtfolder_supp 'mode_1_000000000001_eta_out.nc.trange'], 'LONGITUDE');
for in=1:length(fname)
load(fname{in});
eigvet_map = nan(size(ocean_msk2));
eigvet_map(ocean_msk2) = eigvet(:,1);

hl(in)=plot(lon,eigvet_map,linesty{in});
hold on

%*********table**************************
%%% zonal mean taux
taux = ncread(fnamew{in},'TAUX');
time = ncread(fnamew{in},'TIME');
%detrend
taux_de=nan(size(taux));
for i = 1:341
 for j = 1:61
  reta = squeeze(taux(i,j,:));
  taux_de(i,j,:) = detrend(reta,'linear');
 end
end
%zonal mean
 taux_lont = squeeze(taux_de(:,24,:) + taux_de(:,25,:))/2.;
 snan = isnan(taux_lont);
 taux_lont(snan)=0;
 taux_zmean = sum(taux_lont)./sum(1-snan);
% cor PC1 zmtaux
fname{in}
ny=24;
[r,lags]=xcorr(pc(1,:),taux_zmean,ny,'coeff');
'eigval1,2, cor pc1 and zmtaux at lag=0, max cor'
eigval(1:2)
r(ny+1)
max(r)
lags(r==max(r))

 % significance level
sig = sig_level(taux_zmean,pc(1,:),.95,true);
'sig'
sig
%***********************************

end

line([lon(1) lon(end)],[0 0],'color',[.6 .6 .6])
set(gca,'xtick',lont,'xticklabel',lonstr,'TickDir','out')
set(gca,'ytick',[-9:3:9],'yticklabel',[-9:3:9]);
xlim([130 280])
ylim([-9 9])
hhl=legend(hl(1:4),source{1:4});
set(hhl,'Box','off','location','southeast')
cr=0.1;
title('EOF1 along the equator [cm]')

%%--- Save Figure 9 -------------------------------
cd(rtfolder)
export_fig('Figure9.eps')
export_fig('Figtest.eps')
system(['python ' rtfolder 'mtools/fix_matlab_eps-master/fix_matlab_eps.py Figtest.eps Figure9.eps']);



