%******************************************************************
% This script is for Figure 2
%
% Input:  EP_ustr_anom_mon_mean_NCEP_modelgrid.nc - NCEP 
%	  EP_ERA40_anom_mon_eastward_stress_modelgrid.nc - ERA40
%         eta3d_mm_fit.mat -from two experiments
% Output: Figure2.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/'};

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

fs= 11;
lont = [160 240];
lonstr = {'160\circE'; '120\circW'};

% Postition settings for figure and colorbar
figure('position',[0 0 1000 800])
ax1=axes('position',[.2 .1 .14 .8]);
po1=get(ax1,'position');
leftw = .17; 
cp1=[.105 .1 .025 .8];

%%%--- Data file list -------------------------------------
cd(dtfolder)
fname = {[subfolders{1} 'EP_ustr_anom_mon_mean_NCEP_modelgrid.nc'],...
         [subfolders{2} 'EP_ERA40_anom_mon_eastward_stress_modelgrid.nc'],...
         [subfolders{1} 'eta3d_mm_fit.mat'],...  
         [subfolders{2} 'eta3d_mm_fit.mat']};
tname = {{'(a) NCEP zonal';'wind stress'},...
         {'(b) ERA-40 zonal';'wind stress'},...
         {'(c) NCEP exp.';'sea level'},...
         {'(d) ERA-40 exp.';'sea level'}};    

cp2=[po1(1)+leftw*length(fname),cp1(2:4)]
%%%--- Loop: read, calculate and plot ----------------
lon = ncread([subfolders{1} 'mode_1_000000000001_eta_out.nc.trange'], 'LONGITUDE');

for in=1:length(fname)
 if (in<=2)
  eta = ncread(fname{in},'TAUX');
  time = ncread(fname{in},'TIME');
 else
  load(fname{in});
 end
 eta_lont = squeeze(eta(:,24,:) + eta(:,25,:))/2.;
 
 timenum = datenum('1957-09-01')+time/24./3600.;
 timestr = datestr(timenum,29);
 ttick = timestr;
 for i = 1:length(time)
  ttick(i,:) = [timestr(i,1:8) '01'];
 end  
 t = (datenum(ttick)-datenum('1957-09-01'))*24.*3600;
   
   %%--- plot Figure 2 -------------------------------------   
   if in==1
    axes(ax1);
    contourf(lon,time,eta_lont',[-1.5 -0.06:0.01:0.06 1.5],'edgecolor','no');    
     cmp=cmocean('delta');  % wind cmp  
     colormap(gca,cmp);freezeColors
     caxis([-0.06 0.06])
     
     hct=colorbar('position',cp1);freezeColors
     set(hct,'yAxisLocation','left')
     set(hct,'ylim',[-0.06 0.06],'ytick',[-0.06:0.01:0.06],'TickDir','out')
     ylabel(hct,'Nm^{-2}','fontsize',fs,'fontweight','bold')
        
   elseif in==2
    axes('position',[po1(1)+leftw*(in-1),po1(2:4)]);
    contourf(lon,time,eta_lont',[-1.5 -0.06:0.01:0.06 1.5],'edgecolor','no');
     addpath('/home/xzhu/mtools/cmocean_v1.1/cmocean');
     cmp=cmocean('delta');  % wind cmp     
     colormap(gca,cmp);freezeColors
     caxis([-0.06 0.06])     
   else
    axes('position',[po1(1)+leftw*(in-1),po1(2:4)]);
    contourf(lon,time,eta_lont',[-200 -10:10 200],'edgecolor','no');
     addpath('/home/xzhu/mtools/cmocean_v1.1/cmocean');
     cmp=cmocean('balance');     
     colormap(gca,cmp);freezeColors
     caxis([-10 10])
     
   end  % all contours
   
   if in==length(fname)
     hcs=colorbar('position',cp2);freezeColors
     set(hcs,'yAxisLocation','right')
     set(hcs,'ylim',[-10. 10.],'ytick',[-10.:10.],'TickDir','out')
     ylabel(hcs,'cm','fontsize',fs,'fontweight','bold')
   end
   
   % xyticklabels
   'in= xyticklabels',in
    set(gca,'xtick',lont,'xticklabel',lonstr)
    set(gca,'YTick',t(1:12:end),'YTickLabel',{},'TickDir','out');
    if in==1
    yTick = get(gca,'ytick');
    aa=ttick(1:12:end,1:4);
      for ii=1:length(t(1:12:end))
        if mod(ii,2)==1
        yticklabel = {aa(ii,:)};
        else
        yticklabel = {};
        end      
      text(lon(1)-0.3*(lon(end)-lon(1)),yTick(ii),...
          yticklabel,'VerticalAlignment','bottom')
      end
    end  % if in==1
    title(tname{in},'fontsize',fs)
end

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



