%******************************************************************
% This script is for Figure 8
%
% 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: Figure8.eps
% Note:   After s0_nodetrend_fit_mm_to_aviso_Nmodes_readycoeff.m
%	        s2_detrend_data_EOF_eq_Zhu2017.m
%
% 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/AVISO/'];
dtfolder_supp = [rtfolder 'Greatbatch2018_data/ERA40_exp/'];

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

lont=[120:20:290];
lonstr={'120\circE' '140\circE' '160\circE' '180\circ'...
        '160\circW' '140\circW' '120\circW' '100\circW' '80\circW'};
 
%%--- Read data -------------------------
filename = [dtfolder 'sla3d_EOF_eq.mat']
load(filename);

lat = ncread([dtfolder_supp 'mode_1_000000000001_eta_out.nc.trange'], 'LATITUDE');
lon = ncread([dtfolder_supp 'mode_1_000000000001_eta_out.nc.trange'], 'LONGITUDE');
 
%%--- Loop: plot PC1 and PC2 -------------------
linetxt = {['PC1 ' num2str(eigval(1)) '%'],...
           ['PC2 ' num2str(eigval(2)) '%']};
linesty = {'r' 'b'};
hx = 2;
h1=nan(1,hx);h2=nan(1,hx);

for i = 1:hx

 % 1 Figure 8 subplot 1    
 subplot(2,1,1)
 eigvet_map = nan(size(ocean_msk2));
 eigvet_map(ocean_msk2) = eigvet(:,i);

 h1(i)=plot(lon,eigvet_map,linesty{i});
 line([lon(1) lon(end)],[0 0],'color',[.6 .6 .6])
 hold on

 if i==1
  xlim([130 280])
  ylim([-9 9])
  set(gca,'ytick',[-9:3:9],'yticklabel',[-9:3:9]);
  set(gca,'xtick',lont,'xticklabel',lonstr,'TickDir','out')
 end

 % 2 Figure 8 subplot 2
 subplot(2,1,2)
 h2(i)=plot(time,squeeze(pc(i,:)),linesty{i});
 line([time(1) time(end)],[0 0],'color',[.6 .6 .6])
 hold on

 if i==1
  xlim([time(1) time(end)])
  ylim([-5 5])
  xlabel('year');

  timenum = datenum(time);
  timestr = datestr(timenum,29);
  set(gca,'XTick',time(1:24:end),'XTickLabel',timestr(1:24:end,1:4),'TickDir','out')
  set(gca,'yTick',[-4:2:4],'yticklabel',[-4:2:4],'TickDir','out')
 end

end % end of loop

subplot(2,1,1)
hhl=legend(h1,'EOF1','EOF2')
set(hhl,'Box','off','location','south')
title('(a) EOFs for AVISO along the equator [cm]')

subplot(2,1,2)
hhl=legend(h2,'PC1','PC2')
set(hhl,'Box','off','location','south')
title('(b) normalized principal components time series')

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




