%******************************************************************
% This script is for Figure 6
%
% Input:  eta3d_mm_fit.mat - from four experiments and Zhu et al (2017)
%	  sla3d.mat - AVISO
%	  EP_ustr_anom_mon_mean_NCEP_modelgrid.nc - NCEP 
%	  EP_ERA40_anom_mon_eastward_stress_modelgrid.nc - ERA40
%	  EP_ERA_interim_anom_mon_eastward_stress_modelgrid.nc - ERA-Interim
% Output: Figure6.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 = {'ERA40_exp/' 'NCEP_exp/' 'Zhu2017/' 'AVISO/'};

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

lont=[120:20:280];
lonstr={'120\circE' '140\circE' '160\circE' '180\circ',...
        '160\circW' '140\circW' '120\circW' '100\circW' '80\circW'};

linesty = {'b','b--','k','k--','r-','r--'};

% Setting positions for figure and axes
figure('position',[0 0 550 400])
ax1=axes('position',[.1,.55,.8,.4])
ax2=axes('position',[.1,.05,.8,.4])

%%--- Data file list -------------------------------------
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{3} 'EP_ERA_interim_anom_mon_eastward_stress_modelgrid.nc']};

fnames = {[dtfolder subfolders{1} 'eta3d_mm_fit.mat'],...
	  [dtfolder subfolders{2} 'eta3d_mm_fit.mat'],...
 	  [dtfolder subfolders{3} 'eta3d_mm_fit.mat'],...
	  [dtfolder subfolders{4} 'sla3d.mat']};

%%--- Loop: read data, calculate variance and plot ------------
lon = ncread([dtfolder subfolders{1} 'mode_1_000000000001_eta_out.nc.trange'], 'LONGITUDE');

hlw = nan(1,length(fnamew)+2);
hls = nan(1,length(fnames));
hlw2=hlw;
hls2=hls;
for in=1:length(fnames)+2
 % zonal wind
 if in>length(fnames)
  taux = ncread(fnamew{1},'TAUX');
  taux_lont = squeeze(taux(:,24,:) + taux(:,25,:))/2.;
  varia_taux = nan(341,1);
    for i =1:341
     if (isnan(taux_lont(i,1)))
     else
      if in==length(fnames)+1
       retac = squeeze(taux_lont(i,1:204));       
      else
       retac = squeeze(taux_lont(i,205:end));
      end
       reta = retac - mean(retac);
       varia_taux(i) = var(reta(1:end-1));
     end
    end
 elseif in<=length(fnamew)
  taux = ncread(fnamew{in},'TAUX');
  taux_lont = squeeze(taux(:,24,:) + taux(:,25,:))/2.;
  varia_taux = nan(341,1);
    for i =1:341
     if (isnan(taux_lont(i,1)))
     else
      retac = squeeze(taux_lont(i,:));
      reta = retac - mean(retac);
      varia_taux(i) = var(reta);
     end
    end
 end % if

 % sea level anomalies
 if in<=length(fnames)
  load(fnames{in});
  if in==4
   eta=SLA;
  end
  eta_lont = squeeze(eta(:,24,:) + eta(:,25,:))/2.;
  varia_eta = nan(341,1);mn_eta = nan(341,1); 
    for i =1:341
     if (isnan(eta_lont(i,1)))
     else
     retac = squeeze(eta_lont(i,:));
       reta = retac - mean(retac);
     varia_eta(i) = var(reta);
      mn_eta(i) = mean(reta);
     end
    end
 else
  load(fnames{1});
  eta_lont = squeeze(eta(:,24,:) + eta(:,25,:))/2.;
  varia_eta = nan(341,1);mn_eta = nan(341,1);
    for i =1:341
     if (isnan(eta_lont(i,1)))
     else
      if in==length(fnames)+1
       retac = squeeze(eta_lont(i,1:204));
      else
       retac = squeeze(eta_lont(i,205:end));
      end
       reta = retac - mean(retac);
       varia_eta(i) = var(reta);
       mn_eta(i) = mean(reta);
     end
    end
 end % if

 % plot Figure 6 ---------------------------
 axes(ax1)
 if in==4
 elseif in>length(fnames)
  hlw(in-1)=plot(lon,varia_taux,linesty{in})
  hold on
 else
  in
  varia_taux(100)
  hlw(in)=plot(lon,varia_taux,linesty{in})
  hold on
 end

 if in<=4
  axes(ax2)
  hls(in)=plot(lon,varia_eta,linesty{in})
  hold on
 end

 end % end of loop

axes(ax1)
hhl=legend(hlw,'ERA-40 61-02','NCEP 61-02','ERA-I 93-14','ERA-40 61-77','ERA-40 78-02')
set(hhl,'Box','off')
set(gca,'xtick',lont,'xticklabel',lonstr,'TickDir','out')
xlim([130 280])
ylim([0 6e-4])
title('(a) variance of zonal wind stress along the equator{[N^2}m^{-4}]')

axes(ax2)
hhl=gridLegend(hls,2,'ERA-40 exp. 61-02','NCEP exp. 61-02','ERA-I exp. 93-14','AVISO 93-14',...
           'ERA-40 exp. 61-77','ERA-40 exp. 78-02',...
           'orientation','horizontal','location','north')
set(hhl,'Box','off')
set(gca,'xtick',lont,'xticklabel',lonstr,'TickDir','out')
xlim([130 280])
ylim([0 90]) 
title('(b) variance of sea level along the equator{[cm^2}]')

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


