%******************************************************************
% This script is for Figure 7
%
% Input:  EP_ERA_anom_mon_eastward_stress_19612014_modelgrid.nc - ERA40 (1961-1992) + ERA-Interim (1993-2014)
% Output: Figure7.eps
%
% 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/ERA_combine/'];

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

lont=[120:20:290];
lonstr={'120\circE' '140\circE' '160\circE' '180\circ'...
        '160\circW' '140\circW' '120\circW' '100\circW' '80\circW'};
loct = [91,221]; % 165.5E and 130W

figure('position',[0 0 500 450])

win = 253; % running window 20yr1m = 241m
hlf_win = 126;

%%--- Data file list ------------------------------
fnamew = {[dtfolder 'EP_ERA_anom_mon_eastward_stress_19612014_modelgrid.nc']};

%%--- Read data, calculate running variance and plot ---------------
TIME = ncread(fnamew{1},'TIME');
ln = length(TIME);
lon = ncread(fnamew{1},'LONGITUDE');
lenlon = length(lon);

datestr(TIME)
time=nan(size(TIME));
time = TIME/24/3600+datenum('1957-09-01');
datestr(time)
timestr = datestr(time,29);
  for i = 1:ln
    ttick(i,:) = [timestr(i,1:8) '01'];
  end
xtic = time(hlf_win+1:48:ln-hlf_win);
aa = ttick(hlf_win+1:48:ln-hlf_win,1:4);

linesty={'r','b'};
yl={'left','right'};
hl=nan(1,2);

% 1 loop: read data and calculate var of winds at 165.5E and 130W
for i=1:2
 taux = ncread(fnamew{1},'TAUX',[loct(i),24,1],[1,2,645]);
 taux_max = squeeze(mean(taux,2));

 taux_var=nan(size(taux_max));
 for in=hlf_win+1:ln-hlf_win  
  org = taux_max(in-hlf_win:in+hlf_win);
  org_de = detrend(org,'linear');
  taux_var(in)=var(org_de);  
 end
 
% plot Figure 7 subplot 1
 axes('position',[.1,.57,.8,.34]);hold on
 hl(i)=plot(time(hlf_win+1:ln-hlf_win),taux_var(hlf_win+1:ln-hlf_win),linesty{i})
 set(gca,'yaxislocation',yl{i})
 set(gca,'ycolor',linesty{i})
 xlim([time(hlf_win+1) time(ln-hlf_win)])
 ytic = get(gca,'ytick')
 cr=.1;rt=.2

 set(gca,'box','on');
 set(gca, 'xtick',[],'xticklabel',[])
 set(gca,'color','none');
 eval(['taux_var' num2str(i) '=taux_var;']);

end % end of loop

set(gca,'box','off');
set(gca, 'xtick',xtic,'xticklabel',aa)
xlabel('year')
hhl=legend(hl,'{165.5^o}E, 0^oN','{130^o}W, 0^oN')
set(hhl,'Box','off','location','south')
[r,p]=corrcoef(taux_var1(hlf_win+1:ln-hlf_win),taux_var2(hlf_win+1:ln-hlf_win));
r(1,2)
title({'(a) 21yr running variance of zonal wind stress {[N^2}m^{-4}]'})

% 2 longitude-trend at 165.5E
taux = ncread(fnamew{1},'TAUX',[1,24,1],[341,2,645]);
taux_eq = squeeze(mean(taux,2));

taux_var=nan(size(time));
trd_eq = nan(1,lenlon);
for i=1:lenlon
 if isnan(taux_eq(i,1))
 else
 taux_lon = squeeze(taux_eq(i,:));
 for in=hlf_win+1:ln-hlf_win  
  org = taux_lon(in-hlf_win:in+hlf_win);
  org_de = detrend(org,'linear');
  taux_var(in)=var(org_de);
 end
  org = taux_var(hlf_win+1:ln-hlf_win);  
  de = detrend(org,'linear');
  lin = squeeze(org-de);
  trd_eq(i) = (lin(end)-lin(1))/(ln-1-hlf_win*2)*12; % cm/yr
 end
end

% plot Figure 7 subplot 2
axes('position',[.1,.07,.8,.35],'box','on');hold on
line([lon(1) lon(end)],[0 0],'color',[.6 .6 .6])
hold on
plot(lon,trd_eq,'b')
xlim([140 270])
set(gca,'xtick',lont,'xticklabel',lonstr,'TickDir','out')
title({'(b) trends of 21yr running variance of zonal';
       'wind stress along the equator {[N^2}m^{-4}yr^{-1}]'})

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



    
