clc;clear all;close all;
addpath('/data3/yiwen/matlab_toolbox/seawater/')
pathi='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/'

% file=[pathi 'lead_lag_4_6th_N_bld03.nc']
file=[pathi 'lead_lag_4_6th.nc']
file2=[pathi 'TSRyr_018401-036612.rho_TS_Black_N.volavg.nc']
file3=[pathi 'TSRyr_018401-036612.rho_TS_Black_N.nc']

%-------- Flux data ---------%
n=1:183;
for F1_data=1
load([pathi,'Flux_4_6th_1D.mat']);
DD=DD_M(n);DH=DH_M(n);DT=DT_M(n);
load([pathi,'Heat_Flux_4_6th_1D.mat']);SHF=HF_winter{1}(n);
end
%-------- variables ---------%
PC=ncread(file,'PC');
pcsign= -1; %CTRL==*(-1)
PC1=PC(n,1).*pcsign;

NAO=ncread(file,'NAO').*(-1);
SPG=ncread(file,'SPG');
HBLT=ncread(file,'HBLT');
LABrho=ncread(file,'LABrho');
rhoT=ncread(file,'rhoT');
rhoS=ncread(file,'rhoS');
%----------------------------%
load([pathi,'MLD_4_6th_N.mat'])
MLD=detrend(MLD_Mar,1);

%----------------------------%
file5=[pathi 'BSFyr_018401-036612.nc']
BSF=ncread(file5,'BSF');

%-------- upper ocean rho ----------%
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;
depth=ncread(filename,'z_t');depth_m=depth./100;

file4=[pathi 'TSRyr_018401-036612.nc']
salt=ncread(file4,'SALT');temp=ncread(file4,'TEMP');
PD = sw_dens0(salt,temp);  %kg/m^3

depthIN=depth_m>=0 & depth_m<=150;
index=find(depth_m(depthIN));
depthN=find(depthIN);

area=ncread(filename,'TAREA'); %(cm^2)
loc1=find(LON<=360-37.5 & LON>=360-55 & LAT<=57.5 & LAT>=47.5 );

for k=1:60
areamask{k}=area.*mask(:,:,k);    %2D
areamask_IN{k}=areamask{k}(loc1); %1D
sumarea(k)=squeeze(nansum(areamask_IN{k}));
end

for count=1:size(PD,4)
 for k=depthN(1):depthN(end)
pd_k{count}{k} =PD(:,:,k,count);
pd_IN{count}{k}=pd_k{count}{k}(loc1);
PD_yr(k,count)=nansum(nansum(pd_IN{count}{k}.*areamask_IN{k}))/sumarea(k);
 end
spg_k{count} = BSF(:,:,count);
SPG_IN{count} = spg_k{count}(loc1);
SPG_yr(count)=nansum(nansum(SPG_IN{count}.*areamask_IN{1}))/sumarea(1);
end
SPG_yr=detrend(SPG_yr,1);

%===== 30~200m =====%
PD_yr=PD_yr(6:end,:); 
PD_200m=nanmean(PD_yr,1);
PD_200m=detrend(PD_200m,1);

%----------------------------%
corr1=ncread(file,'corr1');      %NAO,SPG,HBLT,LABrho
corr2=ncread(file,'corr2');      %LABrho,rhot,rhos
regress=ncread(file,'regress3'); %
regress95=ncread(file,'regress95'); 
sigA=ncread(file,'sigcurvesA');
sigB=ncread(file,'sigcurvesB');

rhot_vavg =ncread(file2,'rhot_vavg');
rhos_vavg =ncread(file2,'rhos_vavg');
delr_on_vavg=ncread(file2,'delr_on_vavg');
%=============Fig 8(a)=======================
NUM=length(PC1);
SV_t = zeros(NUM, NUM*2-1);
for shift = -(NUM-1):1:(NUM-1)
    shifted_SV = circshift(PC1, shift);
    SV_t(:, shift + NUM) = shifted_SV;
end
T=SV_t;


Variable={NAO(n) SPG_yr(n) HBLT(n) MLD(n) PD_200m(n) DT SHF} 

for ii=1:size(Variable,2)
 for i=1:size(SV_t,2)
[R,P,RL,RU] = corrcoef(T(:,i),Variable{ii});R1{ii}(i)=R(2);
            
 end
end

%===========================================%
x=repmat(0,1,NUM*2+1);y=linspace(-1,1,NUM*2+1);
lags2=-(NUM-1):1:(NUM-1);
L = length(lags2);
vcrit = sqrt(2)*erfinv(0.95);
lconf = -vcrit/sqrt(L);upconf = vcrit/sqrt(L);

set(0,'defaultfigurecolor','w')
Co=colormap(lines(7));
%=====RHO =====%
for ii=1:size(R1,2)
     L(ii)=plot(lags2,R1{ii},'Linewidth',4,'Color',Co(ii,:));hold on
end
   plot(lags2,[lconf;upconf]*ones(size(lags2)),'r','Linewidth',3,'Linestyle','--');hold on
   plot(x,y,'--','Color',[96 96 96]/225)
   legend([L(1),L(2),L(3),L(4),L(5),L(6),L(7)],...
         {'NAO','SPG','BLD[March]','MLD[March]','RHO_{upper}','D_{in}[thermal]','Surf. Heat Flux'},...
         'Location','southeast','Fontsize',14,'Interpreter','tex');legend boxoff
   ylabel('Corr with PC1');xlabel('lags(years)');
   title('CTRL');
   text(5, .7,'PC1 leads','FontSize',14,'FontWeight','bold')
   set(gca,'Fontsize',18,'Linewidth',2);grid on
   set(gca,'XMinorTick','on','YMinorTick','on','TickDir','out','fontweight', 'bold');ytickformat('%.2f')
   xlim([-15 15]);ylim([-.8 .8]);axis square
