clc ;clear all;close all;
addpath('/data3/yiwen/matlab_toolbox/colorbrewer')

%========== loading data ===========%
pathi='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL'
file1=[pathi 'MOCtot.184-366.eof.atl.nc'] %original data[detrend]

eofmap=ncread(file1,'eof');
lat_depth=ncread(file1,'lat_aux_grid');
depth=ncread(file1,'moc_z');depth_km=depth./1e5;

pcsign=-1;
EOF1=eofmap(:,:,1)*pcsign;
EOF2=eofmap(:,:,2);

s1=EOF1;s2=EOF1;s1(s1<0)=0;s2(s2>0)=0;
s3=EOF2;s4=EOF2;s3(s3<0)=0;s4(s4>0)=0;

lat_loc=lat_depth>=-30;
[xx,yy]=meshgrid(lat_depth(lat_loc),depth_km);

%-------- F1 CTRL --------%
titleStr1=['MOC EOF1= 61.58%'];
titleStr2=['MOC EOF2= 15.45%'];


set(0,'defaultfigurecolor','w')

figure;
pcolor(xx,yy,EOF1(lat_loc,:)');shading interp;hold on
[c1,h1]=contourf(xx,yy,s1(lat_loc,:)',[-1.5:0.2:1.5],'Showtext','on','fill','off','Linestyle','-','LabelSpacing',600);hold on
[c2,h2]=contourf(xx,yy,s2(lat_loc,:)',[-1:0.2:1],'Showtext','on','fill','off','Linestyle','--','LineWidth',2,'LabelSpacing',600);
clabel(c1,h1, 'FontSize', 24, 'FontWeight','bold');clabel(c2, h2, 'FontSize', 24, 'FontWeight','bold');

title(titleStr1);xlabel('latitude');ylabel('depth(km)')
cbar=colorbar('location','eastoutside',...
              'colormap', colormap(flipud(brewermap(30, 'RdYlBu'))))
c=colorbar;c.Title.String = "SV";
set(gca,'YDir','reverse');ytickformat('%.2f')
set(gca,'Fontsize',24,'Linewidth',3,'XMinorTick','on',...
    'YMinorTick','on','TickDir','out','fontweight', 'bold');
caxis([-1.5 1.5]);



figure;
pcolor(xx,yy,EOF2(lat_loc,:)');shading interp;hold on
contourf(xx,yy,s3(lat_loc,:)',[-2:0.1:2],'Showtext','off','fill','off','Linestyle','-');hold on
[c1,h1]=contourf(xx,yy,s3(lat_loc,:)',[-2:0.2:2],'Showtext','on','fill','off','Linestyle','-','LabelSpacing',600);hold on
[c2,h2]=contourf(xx,yy,s4(lat_loc,:)',[-2:0.2:2],'Showtext','on','fill','off','Linestyle','--','LineWidth',2,'LabelSpacing',600);
clabel(c1,h1, 'FontSize', 24, 'FontWeight','bold');clabel(c2, h2, 'FontSize', 24, 'FontWeight','bold');

title(titleStr2);xlabel('latitude');ylabel('depth(km)')
cbar=colorbar('location','eastoutside',...
         'colormap', colormap(flipud(brewermap(30, 'RdYlBu'))))
c=colorbar;c.Title.String = "SV";
set(gca,'YDir','reverse');ytickformat('%.2f')
set(gca,'Fontsize',24,'Linewidth',3,'XMinorTick','on',...
    'YMinorTick','on','TickDir','out','fontweight', 'bold');
caxis([-.8 .8]);


%% ====================================%
% PC timeseries
%====================================%
PC=ncread(file1,'eofts')
PC1=PC(:,1).*(-1);
PC2=PC(:,2);
time=linspace(1958,2018,61);

%=== adding NAO timeseries ====%
pathi='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/'
file=[pathi 'slp.184-366.eof.natl.nc']
NAO=ncread(file,'eofts');
NAO_CTRL=NAO(:,1).*(-1);

color=colormap(lines(10))


subplot(2,1,1)
plot(time,PC1(1:61),'LineStyle','-','Linewidth',3,'Color',color(1,:));hold on
plot(time,PC1(62:122),'LineStyle','-','Linewidth',3,'Color',color(2,:));hold on
plot(time,PC1(123:183),'LineStyle','-','Linewidth',3,'Color',color(3,:));hold on
plot(time,NAO_CTRL(123:183),'-.','Linewidth',3,'Color',color(4,:));hold on
axis([1958 2018 -3 3])
set(gca,'Fontsize',24,'Linewidth',3,'XMinorTick','on',...
    'YMinorTick','on','TickDir','out','fontweight', 'bold');
title('PC1 Time Series');grid on
legend('4th','5th','6th','NAO^+ index[6th]',NumColumns=2);legend boxoff



subplot(2,1,2)
plot(time,PC2(1:61),'Linewidth',3,'Color',color(1,:));hold on
plot(time,PC2(62:122),'Linewidth',3,'Color',color(2,:));hold on
plot(time,PC2(123:183),'Linewidth',3,'Color',color(3,:));hold on
set(gca,'Fontsize',24,'Linewidth',3,'XMinorTick','on',...
    'YMinorTick','on','TickDir','out','fontweight', 'bold');
title('PC2 Time Series');grid on
xlabel('Time(years)','FontSize',20);grid on
legend('4th','5th','6th',NumColumns=1);legend boxoff
axis([1958 2018 -3 3])




