%******************************************************************
% This script is for Figure 5
%
% Input:  eta3d_mm_fit.mat - from four experiments
%	  EP_ustr_anom_mon_mean_NCEP_modelgrid.nc - NCEP 
%	  EP_ERA40_anom_mon_eastward_stress_modelgrid.nc - ERA40
% Output: Figure5.eps
% Note:   After s0_nodetrend_fit_mm_to_aviso_Nmodes_readycoeff.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 = {'NCEP_exp/' 'ERA40_exp/' 'NCEPzonal_exp/' 'ERA40zonal_exp/'};

set(0,'defaultfigurecolor',[1 1 1])
addpath([rtfolder 'mtools/cmocean_v1.1/cmocean']);
addpath([rtfolder 'mtools/freezeColors_v23_cbfreeze']);
addpath([rtfolder 'mtools/altmany-export_fig-9ac0917']);
addpath([rtfolder 'mtools/fix_matlab_eps-master']);

fs= 11;
lont = [150:50:250];
lonstr = {'150\circE';'160\circW';'110\circW'};
latt = [-10:10:10];
latstr = {'10\circS'  'EQ'  '10\circN'};

% setting positions of figure and colorbar
figure('position',[0 0 800 600]);
po1=[.05 .7 .35 .22];
leftw = .42; 
botw = .3;
cp1=[po1(1)+leftw*2-.05, po1(2)-botw+.01, .03, .49];
cp2=[po1(1)+leftw*2-.05, po1(2)-botw*2, .03, .22];

%%--- Data file list ---------------------
fname = {[dtfolder subfolders{1} 'eta3d_mm_fit.mat'],...
         [dtfolder subfolders{2} 'eta3d_mm_fit.mat'],...
         [dtfolder subfolders{3} 'eta3d_mm_fit.mat'],...  
         [dtfolder subfolders{4} 'eta3d_mm_fit.mat'],...
         [dtfolder subfolders{1} 'EP_ustr_anom_mon_mean_NCEP_modelgrid.nc'],...
         [dtfolder subfolders{2} 'EP_ERA40_anom_mon_eastward_stress_modelgrid.nc']};
tname = {{'(a) NCEP exp. sea level trend'},...
         {'(b) ERA40 exp. sea level trend'},...
         {'(c) NCEP-zonal exp. sea level trend'},...
         {'(d) ERA40-zonal exp. sea level trend'},...
         {'(e) NCEP zonal wind stress trend'},...
         {'(f) ERA40 zonal wind stress trend'}};
  
%%--- Read data, calculate trend and plot ---------------
lon = ncread([dtfolder subfolders{1} 'mode_1_000000000001_eta_out.nc.trange'], 'LONGITUDE');
lat = ncread([dtfolder subfolders{1} 'mode_1_000000000001_eta_out.nc.trange'], 'LATITUDE');

for in=1:length(fname)
 % read data
 if (in>=5)
  eta = ncread(fname{in},'TAUX');
 else
  load(fname{in});
 end
 
 trend = nan(341,61);
 lenmon = size(eta,3);
 % trend
 for i = 1:341
  for j = 1:61
   reta = squeeze(eta(i,j,:));
   eta_de(i,j,:) = detrend(reta,'linear');
   lin = squeeze((eta(i,j,:)-eta_de(i,j,:)));
   trend(i,j) = (lin(end)-lin(1))/(lenmon-1)*12; % cm/yr
  end
 end

 % plot 
 axes('position',[po1(1)+leftw*(1-mod(in,2)),po1(2)-botw*(ceil(in/2)-1),po1(3:4)]);
 if (in>=5)
  contourf(lon,lat,trend',[-2e-3 -1e-3:2e-4:1e-3 2e-3],'edgecolor','no');
  cmp=cmocean('delta');  % wind cmp     
  colormap(gca,cmp);freezeColors
  caxis([-1e-3 1e-3])    
  if in==6
   hct=colorbar('position',cp2);freezeColors
   set(hct,'yAxisLocation','right','TickDir','out')
   set(hct,'ylim',[-1e-3 1e-3],'ytick',[-1e-3:5e-4:1e-3])
   xlabel(hct,'Nm^{-2}yr^{-1}','fontsize',fs,'fontweight','bold')
  end
 else 
  contourf(lon,lat,trend',[-1 -0.25:0.05:0.25 -1],'edgecolor','no');
  cmp=cmocean('balance');     %sea level cmp
  colormap(gca,cmp);freezeColors
  caxis([-0.25 0.25])
  if in==4
   hcs=colorbar('position',cp1);freezeColors
   set(hcs,'yAxisLocation','right','TickDir','out')
   set(hcs,'ylim',[-.25 .25],'ytick',[-0.25:0.05:0.25])
   xlabel(hcs,'cm yr^{-1}','fontsize',fs,'fontweight','bold')
  end
 end
 set(gca,'xtick',lont,'xticklabel',lonstr,'TickDir','out')
 set(gca,'ytick',latt,'yticklabel',latstr,'TickDir','out')

 title(tname{in},'fontsize',fs)
end

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





