clc;clear all;close all
addpath('wwcolorbrewer')
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"
pd=ncread(filename,'PD');mask=pd;mask(~isnan(pd))=1;

pathi='....~/G_JRA55_test1/ocn/hist/' % pop2 model archive path


pathi2='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/'

file=[pathi2 'slp.184-366.eof.natl.nc']
file2=[pathi2,'TSRyr_018401-036612.nc']
%====== load PC1 ======%
file1=[pathi2,'lead_lag_4_6th.nc'];
PC=ncread(file1,'PC');
pcsign=-1;PC1=PC(:,1).*pcsign;
MOC_lag2=circshift(PC1,-2);
MOC_lag1=circshift(PC1,-1);

%====== load NAO ======%
NAO=ncread(file,'eof').*(-1);
PC=ncread(file,'eofts');PC1=PC(:,1).*(-1);
PC1_lag1=circshift(PC1,-1);

LON=ncread(file2,'TLONG');LAT=ncread(file2,'TLAT');
temp=ncread(file2,'TEMP');salt=ncread(file2,'SALT');
pd=sw_dens0(salt,temp);
depth=ncread(file2,'z_t');
depth_m=depth./100; %unit=cm

%====== average of upper ocean[200m] temp(C) ======%
SST=squeeze(temp(:,:,1,:));
%====== upper ocean rho =====%
RHO_200=squeeze(nanmean(pd(:,:,50<=depth_m<=150,:),3));


%%====(1)[our study] black box ====%
loc1=LAT>=47.5 & LAT<=57.5 & LON>=360-55 & LON<=360-37.5;
area=ncread(filename,'TAREA');


%===== load SHF ======%
load([pathi2,'Heat_Flux_4_6th.mat'])
ST=1;END=2196; 
month01=ST:12:END;
month02=ST+1:12:END;
month03=ST+2:12:END;
month04=ST+3:12:END;

SHF_shift=circshift(SHF,-1,3);
M={month01 month02 month03 month04}

for ii=1:size(M,2),SHF_month{ii}=SHF_shift(:,:,M{ii});end

for it = 1:size(SHF,3)/12
SHF_yr(:,:,it)=nanmean(SHF(:,:,(it-1)*12+1:(it-1)*12+12),3);
SHF_DJF(:,:,it)=(SHF_month{1}(:,:,it));
end

SHF_lag1=circshift(HF_winter{1},-1);
SHF_lead1=circshift(HF_winter{1},+1);

%===== load MLD ======%
load([pathi2,'MLD_4_6th_N.mat']);clear MLD_DJF
ST=1;END=2196; 

month03=ST+2:12:END;
MLD_Mar=MLD_003(:,:,month03);
%======== correlation =========%
V={SHF_DJF RHO_200 MLD_Mar RHO_200} %variables
lat_real=-80:1:90;
lon_real=0:1:360;
[xq,yq] = meshgrid(lon_real,lat_real);

for ii=1:size(V,2)
 for yr=1:size(SST,3)
V_N{ii}(:,:,yr)=griddata(LON,LAT,V{ii}(:,:,yr),xq,yq); 
 end
V_N2{ii}=circshift(V_N{ii},-180,2);    
V_N2{ii}(:,181:183,:)=[]; 
end

lon2=linspace(-180,180,size(V_N2{1},2));
lat2=linspace(-80,90,size(V_N2{1},1));
%---------- correlation map ------------%
nx=size(V_N2{1},2);ny=size(V_N2{1},1);

mm=cell(1,4);nn=cell(1,4);
for i=1:nx
  for j=1:ny
[R1(j, i), P1(j, i)] = corr(PC1,squeeze(V_N2{1}(j,i,:)));           %NAO vs SHF
[R2(j, i), P2(j, i)] = corr(SHF_lead1',squeeze(V_N2{2}(j,i,:)));    %SHF vs RHO_200
[R3(j, i), P3(j, i)] = corr(HF_winter{1}',squeeze(V_N2{3}(j,i,:))); %SHF vs MLD
[R4(j, i), P4(j, i)] = corr(MOC_lag1,squeeze(V_N2{4}(j,i,:)));      %MOC vs RHO_200
  end
end
P={P1 P2 P3 P4};R={R1 R2 R3 R4}

for ii=1:size(V,2)
[mm{ii} nn{ii}] = find(P{ii}>0.05);
Var=R;Var{ii}(P{ii}<0.05)=nan;
end

[xx,yy]=meshgrid(lon2,lat_real);


Xw=[-55 -55 -37.5 -37.5 -55]
Yw=[47.5 57.5 57.5 47.5 47.5]

set(0,'defaultfigurecolor','w')
time_labels = {'[NAO vs SHF]','[SHF(lead1) vs RHO_{upper}]',...
               '[SHF vs MLD[March]','[MOC(lag1) vs RHO_{upper}]'};

[ha, hfigure]= tight_subplot(2, 2, [.09 .01], [.08 .08], [.08 .08]);

for ii=1:4
  axes(ha(ii));
m_proj('miller','lat',[30,80],'lon',[-80,20]);
m_pcolor(xx,yy,R{ii});shading flat;hold on
m_coast('patch','w','edgecolor','k');
m_grid('xtick',4,'ytick',4,'linestyle','none','tickdir',...
       'out','Linewidth',2,'Fontsize',16);
%---------------------------------------%
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(Xw,Yw,'linewi',2,'color','k');hold on
title({time_labels{ii}},'Interpreter','tex');
if ii==1|ii==3,ylabel('CTRL');end
%---------------------------------------%
colormap(flipud(brewermap(30, 'RdBu')));
caxis([-.7 .7]);
set(gca,'Fontsize',16,'Linewidth',2,'TickDir','out');ytickformat('%.2f')
end

c=colorbar('location','eastoutside','Position', [.88 .09 .01 .85],...
         'colormap', colormap(flipud(brewermap(30,'RdBu'))))
c.Title.String = "correlation";








