%******************************************************************
% This script is for Figure 3 and 4
%
% Input:  psmsl_aviso_sla_[station ID].mat - 1993-2014 tide gague SLA
%	  psmsl_sla_[station ID].mat - 1961-2002 tide gauge SLA 
%	  sla3d_eq_[station ID].nc - 1993-2014 AVISO SLA
%	  eta3d_mm_eq_[station ID].nc - 1961-2002 model-computed SLA
% Output: Figure3.eps
%	  Figure4.eps
% Note:   After s1_compute_psmsl_sla_bothAvisoModel.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/'];
subfolders = {'PSMSL/' 'AVISO/' 'ERA40_exp/' 'NCEP_exp/'};

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

ID_eq = {'1164','1329','1371','1645','927'};
ID_eq = fliplr(ID_eq)  % reverse the order
island_eq = {'RABAUL','KANTON','CHRISTMAS','BALTRA','SANCRIST'};
island_eq = fliplr(island_eq)

lw = 0.8;
legs = 8;
yd = 50;
td = 15;
ytk = [];
ytkl = [];
for ni=1:length(ID_eq)
ytk = [ytk -15+yd*(ni-1):td:15+yd*(ni-1)];
ytkl = [ytkl -15:td:15];
end

hf1=figure('position',[100 100 500 400])
hf2=figure('position',[100 100 600 400])

%%--- Read data and plot --------------------------------------------
for ni=1:length(ID_eq)

 % 1 tide gauge va. AVISO
 fy = 1993.0416;
 ly = 2014.7084;

 load([dtfolder subfolders{1} 'psmsl_aviso_sla_' ID_eq{ni} '.mat'])
 slafilename = [dtfolder subfolders{2} 'sla3d_eq_' ID_eq{ni} '.nc'];
 sla = squeeze(ncread(slafilename,'SLA',[1,1,1],[1,1,length(time)]));

 figure(hf1) 
 hp1 = plot(time,SSL+yd*(ni-1),'k','linewidth',lw);
 line([time(1) time(end)],[yd*(ni-1) yd*(ni-1)],'color','k')
 hold on
 hp2 = plot(time,sla+yd*(ni-1),'color',[.6 .6 .6],'linewidth',lw);
 
 if (ni == 1)
  set(gca, 'xtick',[fy:2:ly],'xticklabel',[floor(fy):2:floor(ly)])
  set(gca, 'ytick',ytk,'yticklabel',ytkl);

  set(gca,'tickdir','out')
  xlim([fy ly])
  ylim([ytk(1)-15 ytk(end)+20])
  xlabel('year')
  ylabel('sea level anomaly [cm]')
 end

 % indicate island name
 yTick=get(gca,'ytick');
 rt = 0.12;
 lll = length(yTick);
 text(time(end)*(1-rt)+time(1)*rt,yTick(floor(lll/5*ni))+yd/5,...
       island_eq{ni},...
          'horizontalAlignment','center','interpreter','tex')

 % 2 tide gauge vs. model
 fy = 1961.0416;
 ly = 2002.6250;

 load([dtfolder subfolders{1} 'psmsl_sla_' ID_eq{ni} '.mat'])
 
 etafilename = {[dtfolder subfolders{3} 'eta3d_mm_eq_' ID_eq{ni} '.nc'],...
                [dtfolder subfolders{4} 'eta3d_mm_eq_' ID_eq{ni} '.nc']};

 figure(hf2)
 hp1 = plot(time,SSL+yd*(ni-1),'k','linewidth',lw);
 line([time(1) time(end)],[yd*(ni-1) yd*(ni-1)],'color','k')
 hold on 

 % loop for NCEP and ERA40 exp.
 hp2 = zeros(1,length(etafilename));
 linesty = {'b' 'r'};
 linlinesty = {'b--' 'r--'};
 r_8597 = nan(1,2);
 
 for nj = 1: length(etafilename)
  eta = squeeze(ncread(etafilename{nj},'ETA',[1,1,1],[1,1,length(time)]));

  hp2(nj) = plot(time,eta+yd*(ni-1),linesty{nj},'linewidth',lw);
    
  eta_de = detrend(eta,'linear');
  lin = eta - eta_de;
  plot(time,lin+yd*(ni-1),linesty{nj},'linewidth',lw/2) % trend
 
 % test print
 island_eq
 ni
 'ERA40  NCEP'
 nj
 (lin(end)-lin(1))/(length(lin)-1)*12  % cm yr-1

  %***** correlations betw. tide gauge and model for 1985-1997
  if ni>1
   fycor = 1985.0416;
   lycor = 1997.9584;
   ify = floor((fycor - fy)*12+1);
   ily = floor((lycor - fy)*12);
   t = time(ify:ily);
   ssl=SSL(ify:ily);
   sssl=ssl(~isnan(ssl));
   stime=t(~isnan(ssl));
   inssl=interp1(stime,sssl,t,'linear','extrap'); 

   SLA4cor(:,6-ni) = inssl;

   [r,p,rl,ru]=corrcoef(inssl,eta(ify:ily));
   r_8597(nj) = r(1,2);  % 0.95 + 

   % test print
%'ssl'
%ssl
%'sssl'
%sssl
'r_8597, rl, ru,p'
[r_8597(nj),rl(1,2),ru(1,2),p(1,2)]

 % significance level
sig = sig_level(inssl,eta(ify:ily),.95,true);
'sig'
sig

  end
  %***********************************************

 end % nj

 if ni==1
 set(gca, 'xtick',[fy:3:ly],'xticklabel',[floor(fy):3:floor(ly)])
 set(gca, 'ytick',ytk,'yticklabel',ytkl);

 set(gca,'tickdir','out')
 xlim([fy ly])
 ylim([ytk(1)-15 ytk(end)+20])
 xlabel('year')
 ylabel('sea level anomaly [cm]') 
 end

 % indicate island name
 yTick=get(gca,'ytick');
 rt = 0.8;
 lll = length(yTick);
 text(time(end)*(1-rt)+time(1)*rt,...
       yTick(floor(lll/5*ni))+yd/5,...
       island_eq{ni},...
          'horizontalAlignment','center','interpreter','tex')
end % ni

%%--- Save Figure 3 and 4 -------------------------------------
cd(rtfolder)
export_fig(hf2,'Figure3.eps')
export_fig(hf2,'Figtest.eps')
system(['python ' rtfolder 'mtools/fix_matlab_eps-master/fix_matlab_eps.py Figtest.eps Figure3.eps']);
export_fig(hf1,'Figure4.eps')
export_fig(hf1,'Figtest.eps')
system(['python ' rtfolder 'mtools/fix_matlab_eps-master/fix_matlab_eps.py Figtest.eps Figure4.eps']);




