%********************************************************************
% This is a pure calculation (no plotting) script: generates PSMSL 
% tide gauge sea level anomalies (relative to monthly mean) over the 
% study period (i) in Zhu et al (2017) for comparison with AVISO 
% (1993-2014), and (ii) in Greatbatch et al (2018) for comparison 
% with model-computed SLA (1961-2002).
%
% Input:  psmsl_aviso_[station ID].mat - 1993-2014 tide gauge sea level
%	  psmsl_[station ID].mat - 1961-2002 tide gauge sea level
% Output: psmsl_aviso_sla_[station ID].mat - sea level anomalies
%	  psmsl_sla_[station ID].mat - sea level anomalies
%
% 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/PSMSL/'];

ID_eq = {'1164','1329','1371','1645','927'};
island_eq = {'RABAUL','KANTON','CHRISTMAS','BALTRA','SANCRIST'};

lw = 0.8;
legs = 8;

%%--- Read sea level and calculate sea level anomalies -----------
cd(dtfolder)
for ni=1:length(ID_eq)

 % 1 tide gauge in aviso time span
 fy = 1993.0416;
 ly = 2014.7084;
 load([dtfolder 'psmsl_aviso_' ID_eq{ni} '.mat']);
   % minus clim mon ssl
  for k =1:12
   sslmon = SSL(k:12:end);  % only one month
   snan = isnan(sslmon);
   sslmon(snan)=0;
   sslmonmean = sum(sslmon)./sum(1-snan);
   sslmona = sslmon - sslmonmean;
   sslmona(snan) = nan;
   SSL(k:12:end) = sslmona;
  end
 %%--- save data ---------------------------------------------
 save(['psmsl_aviso_sla_' ID_eq{ni} '.mat'],'SSL','time')

 % 2 tide gauge in model time span
 fy = 1961.0416;
 ly = 2002.6250;
 load([dtfolder 'psmsl_' ID_eq{ni} '.mat']);
  % minus clim mon ssl
  for k =1:12
   sslmon = SSL(k:12:end);  % only one month
   snan = isnan(sslmon);
   sslmon(snan)=0;
   sslmonmean = sum(sslmon)./sum(1-snan);
   sslmona = sslmon - sslmonmean;
   sslmona(snan) = nan;
   SSL(k:12:end) = sslmona;
  end

 %%--- save data -----------------------------------
 save(['psmsl_sla_' ID_eq{ni} '.mat'],'SSL','time')

end


