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

%========== Model output ===========%
filename="/data3/yiwen/cesm2/JGR_ocean/Model_output/G_JRA55_test1.pop.h.0031-04.nc"
LAT=ncread(filename,'TLAT');LON=ncread(filename,'TLONG');
pd=ncread(filename,'PD');mask=pd;mask(~isnan(pd))=1;

%========== loading data ===========%
pathi='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/'

file1=[pathi1 'lead_lag_4_6th.nc']
pcsign=-1; %CTRL
PC=ncread(file1,'PC');PC1=PC(:,1)*pcsign;

%===== check for NAO & SPG ======%
NAO=ncread(file1,"NAO").*(-1);
SPG=ncread(file1,"SPG");
file5=[pathi 'BSFyr_018401-036612.nc']
BSF=ncread(file5,'BSF');

for lead_lag=1
for i=1:7,PC_lead{i}=circshift(NAO,i,1);end
for i=1:7,PC_lag{i}=circshift(NAO,-i,1);end

nx=320;ny=384;
% ====== !! choose variables ====== %
V1={PC_lead{7} PC_lead{5} PC_lead{3} PC_lead{1}...
    NAO PC_lag{1} PC_lag{3} PC_lag{5} PC_lag{7}}; 

n=size(V1,2);
R1=cell(1,n);P1=cell(1,n);

 for ii=1:n
 for j=1:ny
  for i=1:nx
    [R11(i, j), P11(i,j)] = corr(V1{ii},squeeze(BSF(i,j,:)));  
    R1{ii}=R11;P1{ii}=P11;
  end
 end
end
end

for interpolation=1
%----- interpolation[1X1] -----%
lat_real=-80:1:90;
lon_real=0:1:360;
[xq,yq] = meshgrid(lon_real,lat_real);

for ii=1:n
%(1)
SPG_N{ii}=griddata(LON,LAT,R1{ii},xq,yq);
P1_N{ii}  =griddata(LON,LAT,P1{ii},xq,yq);
%(2)
SPG_N2{ii}=circshift(SPG_N{ii},-180,2);SPG_N2{ii}(:,181:183)=[];
P1_N2{ii}  =circshift(P1_N{ii},-180,2);P1_N2{ii}(:,181:183)=[];
[mm{ii} nn{ii}] = find(P1_N2{ii}>0.05);
V=SPG_N2;V{ii}(P1_N2{ii}>0.05)=nan;
end

end

for region=1
lon2=linspace(-180,180,size(SPG_N2{1},2));
lat2=linspace(-80,90,size(SPG_N2{1},1));
[xx,yy]=meshgrid(lon2,lat2);
%==== study region ====%
Xw2=[-55 -55 -37.5 -37.5 -55]
Yw2=[47.5 57.5 57.5 47.5 47.5]
end


set(0,'defaultfigurecolor','w')
time_labels = {'[NAO lead 7yr]','[NAO lead 5yr]','[NAO lead 3yr]',...
               '[NAO lead 1yr]','[NAO simultaneously]','[NAO lag 1yr]',...
               '[NAO lag 3yrs]','[NAO lag 5yrs]','[NAO lag 7yrs]'};

[ha, hfigure]= tight_subplot(3, 3, [.06 .02], [.15 .05], [.08 .08]);

for ii=1:9
    axes(ha(ii));
m_proj('miller','lat',[30,75],'lon',[-100,30]);
m_pcolor(xx,yy,SPG_N2{ii});shading flat;hold on
m_coast('patch','w','edgecolor','k');
m_grid('linestyle','none','tickdir','out','Linewidth',2,'Fontsize',12);

for i = [1:1:length(mm{ii})],
    m_plot(lon2(nn{ii}(i)),lat_real(mm{ii}(i)),'.',...
          'MarkerEdgeColor','k','MarkerSize',5);hold on
end

m_line(Xw2,Yw2,'linewi',2,'color','k');hold on
colormap(flipud(brewermap(100, 'RdBu')));
title({time_labels{ii}});

set(gca,'Fontsize',10,'Linewidth',2,'TickDir','out');ytickformat('%.2f')
caxis([-.6 .6]) 
end

y=text(-7.5, 1.5,'CTRL: NAO corr. with BSF','FontWeight','bold','FontSize',16)
set(y,'Rotation',90)

c=colorbar('location','southoutside','Position', [0.15 0.05 0.7 0.02],...
         'colormap', colormap(flipud(brewermap(100, 'RdBu'))))
c.Title.String = "correlation";



