% === PURPOSE: REPLOT THE FIGURES FOR MUSCRIPT
% === HISTORY:22023-10-01
% === ATHOUR: EVA
%% LOAD DATA
% % === JRA 
% setup_nctoolbox();
JRADATA = load('DATA\jra55data.mat');
airJRA = JRADATA.air;
hgtJRA = JRADATA.hgt;
fid = 'JRA55/anl_p25.007_hgt.197901_197912';
data = ncdataset(fid);
lon0JRA  = squeeze(data.data('lon'));
latJRA  = squeeze(data.data('lat'));
levelJRA = squeeze(data.data('isobaric'));
% % === NCEP
NCEPDATA = load('DATA/N2andTheta3.mat');
thetaNCEP= NCEPDATA.theta;
N20NCEP = NCEPDATA.N20;
fid = 'NCEPdata/air.mon.mean.nc';
ncdisp(fid)
airNCEP =  ncread(fid,'air');
levelNCEP = ncread(fid,'level');
[xmNCEP,ymNCEP,zmNCEP,tmNCEP] = size(airNCEP);
fid = 'NCEPdata/land.nc';
ncdisp(fid)
landNCEP = ncread(fid,'land');% 1 is land; 0 is ocean
lon0NCEP = double(ncread(fid,'lon'));
latNCEP = double(ncread(fid,'lat'));
lat = latNCEP;
lon0 = lon0NCEP;
% === ERA5
fid = 'ERA5data/era5.ght.mon.mean.nc';
ncdisp(fid)
hgt0 = ncread(fid,'z');% 
hgtERA5 = hgt0./9.8;% 
levelERA5 = double(ncread(fid,'level'));
time = double(ncread(fid,'time'));
latERA5 = ncread(fid,'latitude');
lon0ERA5 = ncread(fid,'longitude');
[xmERA5,ymERA5,zmERA5,tmERA5] = size(hgtERA5);
lon0ERA5(xmERA5/2:xmERA5) =  lon0ERA5(xmERA5/2:xmERA5)-360;
fid = 'ERA5data/era5.temp.mon.mean.nc';
ncdisp(fid)
airERA5 =  ncread(fid,'t');
% ---- POTENTIAL TMEPERATURE
[xmJRA,ymJRA,zmJRA,tmJRA] = size(airJRA);
lon0JRA(xmJRA/2:xmJRA) =  lon0JRA(xmJRA/2:xmJRA)-360;
Rd = 287.058;
gg0 = 9.8 ;% gravity accaleration;
thetaJRA = zeros(xmJRA,ymJRA,zmJRA,tmJRA);% potential temperature;
N20JRA = zeros(xmJRA,ymJRA,zmJRA,tmJRA);% the stability paramenter;
Thro0JRA = zeros(xmJRA,ymJRA,zmJRA,tmJRA);
for k =1:zmJRA
    if k==1
        thetaJRA(:,:,k,:) = airJRA(:,:,k,:).*(1000./levelJRA(k))^0.286;
    elseif k < zmJRA
        thetaJRA(:,:,k+1,:) = airJRA(:,:,k+1,:).*(1000./levelJRA(k+1))^0.286;
        N20JRA(:,:,k,:) = gg0.*(thetaJRA(:,:,k-1,:)-thetaJRA(:,:,k,:))...
            ./(thetaJRA(:,:,k,:).*(hgtJRA(:,:,k-1,:)-hgtJRA(:,:,k,:)));
        Thro0JRA(:,:,k,:) = levelJRA(k)*100/Rd/airJRA(:,:,k,:);
    end
end  
% ==== calculate the N2
Rd = 287.058;
gg0 = 9.8 ;% gravity accaleration;
thetaERA5 = zeros(xmERA5,ymERA5,zmERA5,tmERA5);% potential temperature;
N20ERA5 = zeros(xmERA5,ymERA5,zmERA5,tmERA5);% the stability paramenter;
Thro0ERA5 = zeros(xmERA5,ymERA5,zmERA5,tmERA5);
for k =1:zmERA5
    if k==1
        thetaERA5(:,:,k,:) = airERA5(:,:,k,:).*(1000./levelERA5(k))^0.286;
    elseif k < zmERA5
        thetaERA5(:,:,k+1,:) = airERA5(:,:,k+1,:).*(1000./levelERA5(k+1))^0.286;
        N20ERA5(:,:,k,:) = gg0.*(thetaERA5(:,:,k-1,:)-thetaERA5(:,:,k,:))...
            ./(thetaERA5(:,:,k,:).*(hgtERA5(:,:,k-1,:)-hgtERA5(:,:,k,:)));
        Thro0ERA5(:,:,k,:) = levelERA5(k)*100/Rd/airERA5(:,:,k,:);
    end
end 
clear JRADATA fid Thro0ERA5 Thro0JRA NCEPDATA gg0R
%% === for each month  lower troposphere average omega distribution
mhgtJRA = zeros(xmJRA,ymJRA,zmJRA,12);
dhgtJRA = zeros(xmJRA,ymJRA,zmJRA,tmJRA);
LTN20JRA = squeeze(mean(N20JRA(:,:,27:34,:),3));%700-925
mLTN20JRA = zeros(xmJRA,ymJRA,12);
dLTN20JRA = zeros(xmJRA,ymJRA,tmJRA);
mN20JRA = zeros(xmJRA,ymJRA,zmJRA,12);
dN20JRA = zeros(xmJRA,ymJRA,zmJRA,tmJRA);
mairJRA = zeros(xmJRA,ymJRA,zmJRA,12);
dairJRA = zeros(xmJRA,ymJRA,zmJRA,tmJRA);
mairNCEP = zeros(xmNCEP,ymNCEP,zmNCEP,12);
dairNCEP = zeros(xmNCEP,ymNCEP,zmNCEP,tmNCEP);
mairERA5 = zeros(xmERA5,ymERA5,zmERA5,12);
dairERA5 = zeros(xmERA5,ymERA5,zmERA5,tmERA5);
LTN20ERA5 = squeeze(mean(N20ERA5(:,:,15:16,:),3));%700-850,850-925
mLTN20ERA5 = zeros(xmERA5,ymERA5,12);
dLTN20ERA5 = zeros(xmERA5,ymERA5,tmERA5);
mN20ERA5 = zeros(xmERA5,ymERA5,zmERA5,12);
dN20ERA5 = zeros(xmERA5,ymERA5,zmERA5,tmERA5);
mhgtERA5 = zeros(xmERA5,ymERA5,zmERA5,12);
dhgtERA5 = zeros(xmERA5,ymERA5,zmERA5,tmERA5);
LTN20NCEP = squeeze(mean(N20NCEP(:,:,2:3,:),3));%700-850,850-925
mLTN20NCEP = zeros(xmNCEP,ymNCEP,12);
dLTN20NCEP = zeros(xmNCEP,ymNCEP,tmNCEP);
mN20NCEP = zeros(xmNCEP,ymNCEP,zmNCEP,12);
dN20NCEP = zeros(xmNCEP,ymNCEP,zmNCEP,tmNCEP);
for jt =1:12
    mhgtJRA(:,:,:,jt) = mean(hgtJRA(:,:,:,jt:12:360),4);
    dhgtJRA(:,:,:,jt:12:end) = hgtJRA(:,:,:,jt:12:end)-mhgtJRA(:,:,:,jt);
    mN20JRA(:,:,:,jt) = mean(N20JRA(:,:,:,jt:12:360),4);
    mLTN20JRA(:,:,jt) = mean(LTN20JRA(:,:,jt:12:360),3);
    dN20JRA(:,:,:,jt:12:end) = N20JRA(:,:,:,jt:12:end)-mN20JRA(:,:,:,jt);
    dLTN20JRA(:,:,jt:12:end) = LTN20JRA(:,:,jt:12:end)-mLTN20JRA(:,:,jt);
    mairJRA(:,:,:,jt) = mean(airJRA(:,:,:,jt:12:360),4);
    dairJRA(:,:,:,jt:12:end) = airJRA(:,:,:,jt:12:end)-mairJRA(:,:,:,jt);
    mairNCEP(:,:,:,jt) = mean(airNCEP(:,:,:,jt:12:360),4);
    dairNCEP(:,:,:,jt:12:end) = airNCEP(:,:,:,jt:12:end)-mairNCEP(:,:,:,jt);
    mairERA5(:,:,:,jt) = mean(airERA5(:,:,:,jt:12:360),4);
    dairERA5(:,:,:,jt:12:end) = airERA5(:,:,:,jt:12:end)-mairERA5(:,:,:,jt);
    mhgtERA5(:,:,:,jt) = mean(hgtERA5(:,:,:,jt:12:360),4);
    dhgtERA5(:,:,:,jt:12:end) = hgtERA5(:,:,:,jt:12:end)-mhgtERA5(:,:,:,jt);
    mN20ERA5(:,:,:,jt) = mean(N20ERA5(:,:,:,jt:12:360),4);
    mLTN20ERA5(:,:,jt) = mean(LTN20ERA5(:,:,jt:12:360),3);
    dN20ERA5(:,:,:,jt:12:end) = N20ERA5(:,:,:,jt:12:end)-mN20ERA5(:,:,:,jt);
    dLTN20ERA5(:,:,jt:12:end) = LTN20ERA5(:,:,jt:12:end)-mLTN20ERA5(:,:,jt);
    mN20NCEP(:,:,:,jt) = mean(N20NCEP(:,:,:,jt:12:360),4);
    mLTN20NCEP(:,:,jt) = mean(LTN20NCEP(:,:,jt:12:360),3);
    dN20NCEP(:,:,:,jt:12:end) = N20NCEP(:,:,:,jt:12:end)-mN20NCEP(:,:,:,jt);
    dLTN20NCEP(:,:,jt:12:end) = LTN20NCEP(:,:,jt:12:end)-mLTN20NCEP(:,:,jt);
end
T2014 = datestr((datenum('19000101','yyyymmdd')+time./24),'yyyymm');

% === winter multiyear average lower atmosphere stability
% === the first is the 1980 winter 2023 winter
data0 = squeeze(LTN20JRA(:,:,12:12:tmJRA-1));
data1 = squeeze(LTN20JRA(:,:,13:12:tmJRA-1));
data2 = squeeze(LTN20JRA(:,:,14:12:tmJRA-1));
data3 = squeeze(LTN20JRA(:,:,15:12:tmJRA-1));
WinLtN20JRA = (data0+data1+data2+data3)./4;% average monthly
mWintrLTN20JRA = mean(WinLtN20JRA,3);
%
data0 = squeeze(LTN20ERA5(:,:,12:12:tmERA5-1));
data1 = squeeze(LTN20ERA5(:,:,13:12:tmERA5-1));
data2 = squeeze(LTN20ERA5(:,:,14:12:tmERA5-1));
data3 = squeeze(LTN20ERA5(:,:,15:12:tmERA5-1));
WinLtN20ERA5 = (data0+data1+data2+data3)./4;% average monthly
mWintrLTN20ERA5 = mean(WinLtN20ERA5,3);
%
data0 = squeeze(LTN20NCEP(:,:,12:12:tmNCEP-1));
data1 = squeeze(LTN20NCEP(:,:,13:12:tmNCEP-1));
data2 = squeeze(LTN20NCEP(:,:,14:12:tmNCEP-1));
data3 = squeeze(LTN20NCEP(:,:,15:12:tmNCEP-1));
WinLtN20NCEP = (data0+data1+data2+data3)./4;% average monthly
mWintrLTN20NCEP = mean(WinLtN20NCEP,3);
% === summer multiyear average lower atmosphere stability
% === the first is the 1980 summer 2021summer
data0 = squeeze(LTN20JRA(:,:,6+12:12:tmJRA-2));
data1 = squeeze(LTN20JRA(:,:,7+12:12:tmJRA-2));
data2 = squeeze(LTN20JRA(:,:,8+12:12:tmJRA-2));
SummerLtN20JRA = (data0+data1+data2)./3;% average monthly
mSummertrLTN20JRA = mean(SummerLtN20JRA,3);
%
data0 = squeeze(LTN20ERA5(:,:,6+12:12:tmERA5-2));
data1 = squeeze(LTN20ERA5(:,:,7+12:12:tmERA5-2));
data2 = squeeze(LTN20ERA5(:,:,8+12:12:tmERA5-2));
SummerLtN20ERA5 = (data0+data1+data2)./3;% average monthly
mSummertrLTN20ERA5 = mean(SummerLtN20ERA5,3);
%
data0 = squeeze(LTN20NCEP(:,:,6+12:12:tmNCEP-2));
data1 = squeeze(LTN20NCEP(:,:,7+12:12:tmNCEP-2));
data2 = squeeze(LTN20NCEP(:,:,8+12:12:tmNCEP-2));
SummerLtN20NCEP = (data0+data1+data2)./3;% average monthly
mSummertrLTN20NCEP = mean(SummerLtN20NCEP,3);

% === winter trend between 1979 and 2021
trendAllWinltN2JRA = ones(xmJRA,ymJRA).*nan;
r_matriAllWinltN2JRA = ones(xmJRA,ymJRA)*nan;
errWinltN2JRA = ones(xmJRA,ymJRA)*nan;
trendAllWinltN2ERA5 = ones(xmERA5,ymERA5).*nan;
r_matriAllWinltN2ERA5 = ones(xmERA5,ymERA5)*nan;
errWinltN2ERA5 = ones(xmERA5,ymERA5)*nan;
trendAllWinltN2NCEP = ones(xmNCEP,ymNCEP).*nan;
r_matriAllWinltN2NCEP = ones(xmNCEP,ymNCEP)*nan;
errWinltN2NCEP = ones(xmNCEP,ymNCEP)*nan;
for i = 1:xmJRA
    for j = 1:25
        % === the each point time series: each winter between 1979-2020
        % === 1979-2021 JRA55
        ln2m = squeeze(WinLtN20JRA(i,j,:));
        lmm = length(ln2m);
        % -- trend
        xx = (1:lmm)';
        aa= polyfit(xx,ln2m,1);
        yy = aa(1).*xx + aa(2);
        trendAllWinltN2JRA(i,j) = aa(1);
        % P TEST
        errWinltN2JRA(i,j) = (sum((ln2m - yy).*(ln2m - yy))./42)^0.5;
        r_matriAllWinltN2JRA(i,j)  = corr(xx,ln2m);
        % === 1979-2021 ERA5
        ln2m = squeeze(WinLtN20ERA5(i,j,:));
        lmm = length(ln2m);
        % -- trend
        xx = (1:lmm)';
        aa= polyfit(xx,ln2m,1);
        yy = aa(1).*xx + aa(2);
        trendAllWinltN2ERA5(i,j) = aa(1);
        % P TEST
        errWinltN2ERA5(i,j) = (sum((ln2m - yy).*(ln2m - yy))./42)^0.5;
        r_matriAllWinltN2ERA5(i,j)  = corr(xx,ln2m);
        ln2m = squeeze(WinLtN20NCEP(i,j,:));
        % === 1979-2021
        lmm = length(ln2m);
        % -- trend
        xx = (1:lmm)';
        aa= polyfit(xx,ln2m,1);
        yy = aa(1).*xx + aa(2);
        trendAllWinltN2NCEP(i,j) = aa(1);
        % 
        errWinltN2NCEP(i,j) = (sum((ln2m - yy).*(ln2m - yy))./42)^0.5;
        r_matriAllWinltN2NCEP(i,j)  = corr(xx,ln2m);
    end
end

% === summer N2 trend between 1979 and 2021
trendAllSummerltN2JRA = ones(xmJRA,ymJRA).*nan;
r_matriAllSummerltN2JRA = ones(xmJRA,ymJRA)*nan;
errSumltN2JRA = ones(xmJRA,ymJRA)*nan;
trendAllSummerltN2ERA5 = ones(xmERA5,ymERA5).*nan;
r_matriAllSummerltN2ERA5 = ones(xmERA5,ymERA5)*nan;
errSumltN2ERA5 = ones(xmERA5,ymERA5)*nan;
trendAllSummerltN2NCEP = ones(xmNCEP,ymNCEP).*nan;
r_matriAllSummerltN2NCEP = ones(xmNCEP,ymNCEP)*nan;
errSumltN2NCEP = ones(xmNCEP,ymNCEP)*nan;
for i = 1:xmJRA
    for j = 1:25
        % === the each point time series: each Summer between 1979-2020
        % === 1979-2021 JRA55
        ln2m = squeeze(SummerLtN20JRA(i,j,:));
        lmm = length(ln2m);
        % -- trend
        xx = (1:lmm)';
        aa= polyfit(xx,ln2m,1);
        yy = aa(1).*xx + aa(2);
        trendAllSummerltN2JRA(i,j) = aa(1);
        errSumltN2JRA(i,j) = (sum((ln2m - yy).*(ln2m - yy))./42)^0.5;
        r_matriAllSummerltN2JRA(i,j)  = corr(xx,ln2m);
        % === 1979-2021 ERA5
        ln2m = squeeze(SummerLtN20ERA5(i,j,:));
        lmm = length(ln2m);
        % -- trend
        xx = (1:lmm)';
        aa= polyfit(xx,ln2m,1);
        yy = aa(1).*xx + aa(2);
        trendAllSummerltN2ERA5(i,j) = aa(1);
        errSumltN2ERA5(i,j) = (sum((ln2m - yy).*(ln2m - yy))./42)^0.5;
        r_matriAllSummerltN2ERA5(i,j)  = corr(xx,ln2m);
        % === 1979-2021 NCEP
        ln2m = squeeze(SummerLtN20NCEP(i,j,:));
        lmm = length(ln2m);
        % -- trend
        xx = (1:lmm)';
        aa= polyfit(xx,ln2m,1);
        yy = aa(1).*xx + aa(2);
        trendAllSummerltN2NCEP(i,j) = aa(1);
        errSumltN2NCEP(i,j) = (sum((ln2m - yy).*(ln2m - yy))./42)^0.5;
        r_matriAllSummerltN2NCEP(i,j)  = corr(xx,ln2m);
   end
end
diffLtN2SumJRA = trendAllSummerltN2JRA.*(2021-1979);
diffLtN2WinJRA = trendAllWinltN2JRA.*(2021-1979);
diffLtN2SumERA5 = trendAllSummerltN2ERA5.*(2021-1979);
diffLtN2WinERA5 = trendAllWinltN2ERA5.*(2021-1979);
diffLtN2SumNCEP = trendAllSummerltN2NCEP.*(2021-1979);
diffLtN2WinNCEP = trendAllWinltN2NCEP.*(2021-1979);
%% mask
maskNor0 = mWintrLTN20NCEP.*nan;
lon0NCEP(xmNCEP/2:xmNCEP) =  lon0NCEP(xmNCEP/2:xmNCEP)-360;
id1 = find(lon0NCEP==160);
id2 = find(lon0NCEP==-125);

maskNor0(1:id1,:) = 1;
maskNor0(id2:xmNCEP,:) = 5;
maskNor0(:,24:ymNCEP) = nan;
maskNor0(:,1:14) = nan;
%
maskNor = mWintrLTN20NCEP.*nan;
maskNor(mWintrLTN20NCEP<=0.000134) = 2;
maskNor(109:130,7:17) = nan;
maskNor(:,22:ymNCEP) = nan;
[Plg0,~]=meshgrid(lon0NCEP,latNCEP);
Plg0(maskNor'==2) = 2000;
Plg0(:,35:109) = nan;
Plg0(18:ymNCEP,:) = nan;
Plg0(Plg0<400) = nan;
maskNor(Plg0'==2000) = 1;% for southern iceland sea
[Plg0,Plt0]=meshgrid(lon0NCEP,latNCEP);
Plg0(maskNor'==1) = 2000;
Plg0(:,35:114) = nan;
Plg0(18:ymNCEP,:) = nan;
Plg0(Plt0<64) = nan;
maskNor(Plg0'==2000) = 3;%for norweigain sea
maskNor(:,1:5) = 4;% north of 80N
maskNor(landNCEP==1) = nan;
maskNor(maskNor0==6) = 6;% region-c
%% === fig1 === plot the multiyear average for each month
figure;
set(gcf,'color','white','Units','centimeters','Position',[3 3 18 22.5]);
m_proj('stereographic','lat',90,'long',-10,'radius',35);
Nh=4; Nw=3;
% ha = tight_subplot(Nh, Nw, gap, marg_h, marg_w);%
ha = tight_subplot(Nh, Nw, [0.03,0.04], [0.09,0.04], [0.02,0.02]);%
f_id = 1;
lon = cat(1,lon0NCEP,lon0NCEP(1,:));
[Plg,Plt]=meshgrid(lon,latNCEP);
lat = latNCEP;
MonthStr = ['a JAN';'b FEB';'c MAR';'d APR';'e MAY';'f JUN';'g JUL';...
    'h AUG';'i SEP';'j OCT';'k NOV';'l DEC'];
for tt = 1:12
    %     clf
%     data0 = squeeze(mLTN20JRA(:,:,tt));
    data0 = squeeze(mLTN20NCEP(:,:,tt));
%     data0 = squeeze(mLTN20ERA5(:,:,tt));
    data = cat(1,data0,data0(1,:));
    axes(ha(tt));
    m_contourf(Plg,Plt,data',30,'linestyle','none');shading flat;%colormap(map);
    hold on
    %     hold on;m_contour(Plg,Plt,data',[mid0 mid0],'linestyle','-','color','white');
    %     hold on;m_contour(Plg,Plt,data',[mid1 mid1],'linestyle','-','color','black');
    shading flat;colormap(m_colmap('jet'));
%     m_coast('color','k','linewidth',1.5);
    m_coast('patch',[.7 .7 .7],'edgecolor',[.7 .7 .7]);
    m_grid('box','fancy','FontSize',8,'FontWeight','bold','xtick',-180:40:360,...
        'ytick',30:10:90,'linest',':','tickdir','in');% add a standard colorbar.
    text(-0.56,0.65,{MonthStr(tt,:),' '},'fontangle','normal','fontweight','bold','FontSize',8);
    caxis([1 4]*0.0001);
    dataP = cat(1,maskNor,maskNor(1,:));
    [x,y] =  find(dataP==6);
    kk = boundary(x,y);
    [x,y] =  find(dataP==1);
    kk = boundary(x,y);
    m_plot(lon(x(kk)),lat(y(kk)),'y','linewidth',1)
    [x,y] =  find(dataP==4);
    kk = boundary(x,y);
    m_plot(lon(x(kk(149:270))),lat(y(kk(149:270))),'Color',[0.6 0.2 1],'linewidth',1)
    dataP2 = circshift(dataP,70,1);
    lon2 = circshift(lon,70,1);
    [x,y] =  find(dataP2==3);
    kk = boundary(x,y);
    m_plot(lon2(x(kk)),lat(y(kk)),'r','linewidth',1)
%     lon2 = [0:360];
%     lat2 = lon2.*0+80;
%     m_plot(lon2,lat2,'y','linewidth',1)
    %
    if tt==12
        h2=colorbar('h');%horizontal
        set(h2,'position',[0.06 0.05 0.88 0.02]); % ??position:?left??bottom??width??height
        set(get(h2,'label'),'string','s^{-2}','FontWeight','bold',...
            'rotation',0,'position',[0.00041,0.1])
        h1=axes('position',[0.05 0.05 0.9 0.9]);
%         title(h1,['Fig. 1 Climatological monthly static stability (N^2) ' ...
%             'of the Arctic lower atmosphere averaged over 1979–2021 (unit: s^{-2}) based on JRA55'],'fontangle','normal','FontSize',8)
%         title(h1,['Fig. 1 Climatological monthly static stability (N^2) ' ...
%             'of the Arctic lower atmosphere averaged over 1979–2021 (unit: s^{-2}) based on NCEP'],'fontangle','normal','FontSize',8)
%         title(h1,['Fig. 1 Climatological monthly static stability (N^2) ' ...
%             'of the Arctic lower atmosphere averaged over 1979–2021 (unit: s^{-2}) based on ERA5'],'fontangle','normal','FontSize',8)
        axis off
    end
end
set(gca,'FontSize',8,'FontWeight','bold','LineWidth',1.5,'FontName','Helvetica');
% sav_name=['npj/final submission/JRA55FIG1withoutLand'];
sav_name=['npj/final submission/NCEPFIG1withoutLand'];
% sav_name=['npj/final submission/ERA5FIG1withoutLand'];
print ('-r600','-depsc2',(sav_name))
print ('-r600','-dpdf',(sav_name))
%% === new fig2 for N2 and 1000hpa air temperature
% without and with running average
aoi0 = textread('NCEPdata/aoi.txt');
f1 = find(aoi0(:,1)==1979,1);
aoi = aoi0(f1:f1+tmERA5-1,:,:);
% === the origional time seires AOI
AOIor = aoi(1:tmERA5,3);
% --- running average
nd = 6;
nd2 = 13;
AOIorRA0 = running_average(AOIor,nd,nd,nd2);
AOIorRA = AOIorRA0(nd+1:tmERA5-nd);
txx = aoi(nd+1:tmERA5-nd,1)+aoi(nd+1:tmERA5-nd,2)./12;
% === the cc between aoi and N2
% --- running average
nd = 6;
nd2 = 13;
AOIorRA0 = running_average(AOIor,nd,nd,nd2);
AOIorRA = AOIorRA0(nd+1:tmERA5-nd);
txx = aoi(nd+1:tmERA5-nd,1)+aoi(nd+1:tmERA5-nd,2)./12;
% === the cc between aoi and N2
ccN2TsurfRANCEP = zeros(xmNCEP,ymNCEP);
CCN2TsurfNCEP = zeros(xmNCEP,ymNCEP);
ccN2TsurfRAJRA = zeros(xmJRA,ymJRA);
CCN2TsurfJRA = zeros(xmJRA,ymJRA);
ccN2TsurfRAERA5 = zeros(xmERA5,ymERA5);
CCN2TsurfERA5 = zeros(xmERA5,ymERA5);
for i = 1:xmNCEP
    for j = 1:ymNCEP
        %NCEP
        ln1mAll =  (squeeze(dLTN20NCEP(i,j,1:tmERA5)));
        ln1mAllRA0 = running_average(ln1mAll,nd,nd,nd2);
        ln1mAllRA = ln1mAllRA0(nd+1:tmERA5-nd);
        %
        ln2mAll =  (squeeze(dairNCEP(i,j,1,1:tmERA5)));
        ln2mAllRA0 = running_average(ln2mAll,nd,nd,nd2);
        ln2mAllRA = ln2mAllRA0(nd+1:tmERA5-nd);
        CCN2TsurfNCEP(i,j) = corr(ln1mAll,ln2mAll);
        ccN2TsurfRANCEP(i,j) = corr(ln1mAllRA',ln2mAllRA');
        % JRA
        ln1mAll =  (squeeze(dLTN20JRA(i,j,1:tmERA5)));
        ln1mAllRA0 = running_average(ln1mAll,nd,nd,nd2);
        ln1mAllRA = ln1mAllRA0(nd+1:tmERA5-nd);
        %
        ln2mAll =  (squeeze(dairJRA(i,j,37,1:tmERA5)));
        ln2mAllRA0 = running_average(ln2mAll,nd,nd,nd2);
        ln2mAllRA = ln2mAllRA0(nd+1:tmERA5-nd);
        CCN2TsurfJRA(i,j) = corr(ln1mAll,ln2mAll);
        ccN2TsurfRAJRA(i,j) = corr(ln1mAllRA',ln2mAllRA');
        % ERA5
        ln1mAll =  (squeeze(dLTN20ERA5(i,j,1:tmERA5)));
        ln1mAllRA0 = running_average(ln1mAll,nd,nd,nd2);
        ln1mAllRA = ln1mAllRA0(nd+1:tmERA5-nd);
        %
        ln2mAll =  (squeeze(dairERA5(i,j,17,1:tmERA5)));
        ln2mAllRA0 = running_average(ln2mAll,nd,nd,nd2);
        ln2mAllRA = ln2mAllRA0(nd+1:tmERA5-nd);
        CCN2TsurfERA5(i,j) = corr(ln1mAll,ln2mAll);
        ccN2TsurfRAERA5(i,j) = corr(ln1mAllRA',ln2mAllRA');

    end
end
%% ====== FIG3a ====
figure;
set(gcf,'color','white','Units','centimeters','Position',[3 3 6 6.5]);
Nh=1; Nw=1; gap=0.04; 
ha = tight_subplot(Nh, Nw, [0.04,0.02], [0.15,0.08], [0.1,0.1]);%
m_proj('stereographic','lat',90,'long',-10,'radius',35);
lon = cat(1,lon0NCEP,lon0NCEP(1,:));
[Plg,Plt]=meshgrid(lon,latNCEP);
axes(ha(1));%
data0 = ccN2TsurfRAERA5;
data = cat(1,data0,data0(1,:));
m_contourf(Plg,Plt,data',20,'linestyle','none');shading flat;%colormahold on
shading flat;colormap(m_colmap('divergent'));hold on
hold on;m_contour(Plg,Plt,data',[0 0],'linestyle','--','color','black');
hold on;m_contour(Plg,Plt,data',[-0.4227 -0.4227],'linestyle','-','color','black');
% m_coast('color','k','linewidth',1.5);
m_coast('patch',[.7 .7 .7],'edgecolor',[.7 .7 .7]);
m_grid('box','fancy','FontSize',8,'FontWeight','bold','xtick',[-180:40:360],...
                'ytick',[30:10:90],'linest',':','tickdir','in');% add a standard colorbar.
caxis([-1 1])
% ===========================================
dataP = cat(1,maskNor,maskNor(1,:));
[x,y] =  find(dataP==6);
kk = boundary(x,y);
[x,y] =  find(dataP==1);
kk = boundary(x,y);
m_plot(lon(x(kk)),lat(y(kk)),'y','linewidth',1)
[x,y] =  find(dataP==4);
kk = boundary(x,y);
m_plot(lon(x(kk(149:270))),lat(y(kk(149:270))),'Color',[0.6 0.2 1],'linewidth',1)
dataP2 = circshift(dataP,70,1);
lon2 = circshift(lon,70,1);
[x,y] =  find(dataP2==3);
kk = boundary(x,y);
m_plot(lon2(x(kk)),lat(y(kk)),'r','linewidth',1)
% ===========================================
% a
text(-0.6,0.65 ,{'a'},...
    'LineStyle','none',...
    'FontWeight','bold',...
    'FontSize',8,...
    'FontName','Helvetica');
%CC of N^2 and surface air temperature
h2=colorbar('h');%horizontal
set(h2,'position',[0.1 0.07 0.8 0.04]); % ??position:?left??bottom??width??height
% title({'Summer'},'FontSize',8,'position',[0,0.67]);
sav_name=['npj/final submission/Fig3aCCN2T1000ERA5'];
print ('-r600','-depsc2',(sav_name))
print ('-r600','-dpdf',(sav_name))
print ('-r600','-djpeg',(sav_name))

%%
lmordN2PolNCEP = ones(xmJRA*7,tmERA5).*nan;
lmordair1qPolNCEP = ones(xmJRA*7,tmERA5).*nan;

lmordN2PolERA5 = ones(xmJRA*7,tmERA5).*nan;
lmordair1qPolERA5 = ones(xmJRA*7,tmERA5).*nan;

lmordN2PolJRA = ones(xmJRA*7,tmERA5).*nan;
lmordair1qPolJRA = ones(xmJRA*7,tmERA5).*nan;

kk3 = 1;
for i = 1:xmJRA
    for j = 1:25
        % === calculate the average amplitude without the trend
        ln1mNCEP = squeeze(dLTN20NCEP(i,j,1:tmERA5));
        ln2mNCEP = squeeze(dairNCEP(i,j,1,1:tmERA5));
        ln1mERA5 = squeeze(dLTN20ERA5(i,j,1:tmERA5));
        ln2mERA5 = squeeze(dairERA5(i,j,17,1:tmERA5));
        ln1mJRA = squeeze(dLTN20JRA(i,j,1:tmERA5));
        ln2mJRA = squeeze(dairJRA(i,j,37,1:tmERA5));
        if (maskNor(i,j)==4)
            lmordN2PolNCEP(kk3,:) = ln1mNCEP;
            lmordair1qPolNCEP(kk3,:) = ln2mNCEP;
            lmordN2PolERA5(kk3,:) = ln1mERA5;
            lmordair1qPolERA5(kk3,:) = ln2mERA5;
            lmordN2PolJRA(kk3,:) = ln1mJRA;
            lmordair1qPolJRA(kk3,:) = ln2mJRA;
            kk3 = kk3+1;
        end
    end
end
%% === fig 3 bc
txx = aoi(nd+1:tmERA5-nd,1)+aoi(nd+1:tmERA5-nd,2)./12;
figure;
set(gcf,'color','white','Units','centimeters','Position',[3 3 6 6.5]);
Nh=2; Nw=1; gap=0.04; 
ha = tight_subplot(Nh, Nw, [0.04,0.01], [0.17,0.06], [0.18,0.13]);%
axes(ha(1));
lmAgdN2PolNCEP = nanmean(lmordN2PolNCEP);
lmAgdair1qPolNCEP = nanmean(lmordair1qPolNCEP);
lm1mid = running_average(lmAgdN2PolNCEP,nd,nd,nd2); 
lm2mid = running_average(lmAgdair1qPolNCEP,nd,nd,nd2); 
lm1 = lm1mid(nd+1:tmERA5-nd);
lm2 = lm2mid(nd+1:tmERA5-nd);
[hAx,hLine1,hLine2] = plotyy(txx,lm1,txx,lm2);
hLine1.LineWidth = 1.5;
hLine2.LineWidth = 1.5;
xlim(hAx(1),[1980,2021])
xlim(hAx(2),[1980,2021])
% ylim(hAx(1),[-0.00002 0.00002])
% ylim(hAx(2),[-0.02,0.02])
ylabel(hAx(1),{'N^2 anomaly (s^{-2})'})
ylabel(hAx(2),'SAT anomaly (K)')
set(hAx(1),'FontSize',8,'LineWidth',1,'TickDir',...
    'in','XTick',[1980:5:2021],'XMinorTick','on','XTickLabel',...
    [],'TickLength',[0.025 0.05],...
      'YMinorTick','on','FontName','Helvetica');
set(hAx(2),'FontSize',8,'LineWidth',1,'TickDir',...
    'in','XTick',[1980:5:2021],'XMinorTick','on','TickLength',[0.025 0.05],'XTickLabel',...
    [],...
    'YMinorTick','on');text('Parent',gca,'VerticalAlignment','baseline','FontWeight','bold',...
    'FontSize',8,'FontName','Helvetica',...
    'String',{'Polar Region'},...
    'Position',[1980.7207729383 3e-05 0]);
text('Parent',gca,'VerticalAlignment','baseline','FontWeight','bold',...
    'FontSize',8,'FontName','Helvetica',...
    'String',{'b'},...
    'Position',[1981.20891518602 -9*0.00001 0]);
hold on;plot(txx,txx*0,'k-')
%lm1,txx,lm2
axes(ha(2));
nd = 6;
tmm = length(lm2);
syCC2 = zeros(tmm,1);
a1 = mean(lm1(nd:tmm-nd));
a2 = mean(lm2(nd:tmm-nd));
for kk = nd+1:tmm-nd-1
    xk = lm1(kk-nd:kk+nd);
    yk = lm2(kk-nd:kk+nd);
%     a1 = mean(xk);
%     a2 = mean(yk);
    dxk = xk-a1;
    dyk = yk-a2;
    ntr = sum(dxk.*dyk);
    dtr = (sum(dxk.*dxk).*sum(dyk.*dyk))^0.5;
    syCC2(kk) = ntr./dtr;
end
area(txx,syCC2)
hold on;plot(txx,txx*0+0.35,'k--')
hold on;plot(txx,txx*0-0.35,'k--')
set(gca,'FontSize',8,'LineWidth',1,'FontName','Helvetica');
set(gca,'FontName','Helvetica','FontSize',8,'TickDir','in',...
    'XMinorTick','on','YMinorTick','on','TickLength',[0.025 0.05],'XTick',[1980:5:2021],'XTickLabel',...
    {'1980','1985','1990','1995','2000','2005','2010','2015','2020'});
ylabel({'Correlation ','coefficient'})
xlabel('Year')
xlim([1980,2021])
text('Parent',gca,'VerticalAlignment','baseline','FontWeight','bold',...
    'FontSize',8,'FontName','Helvetica',...
    'String',{'c'},...
    'Position',[1980.83601087494 -0.867539046010974 0]);
sav_name=['npj/final submission/Fig3bcRCCN2T1000NCEP'];
print ('-r600','-depsc2',(sav_name))
print ('-r600','-dpdf',(sav_name))
print ('-r600','-djpeg',(sav_name))
%% ==== fig 2 the trend of N2 between 1979-2020
% === change the unit with s-2
figure;
set(gcf,'color','white','Units','centimeters','Position',[3 3 12 16]);
Nh=3; Nw=2; gap=0.04; 
ha = tight_subplot(Nh, Nw, [0.04,0.02], [0.06,0.06], [0.02,0.1]);%
m_proj('stereographic','lat',90,'long',-10,'radius',35);
lon = cat(1,lon0NCEP,lon0NCEP(1,:));
[Plg,Plt]=meshgrid(lon,latNCEP);
%
axes(ha(1));
data0 = diffLtN2SumNCEP;
data = cat(1,data0,data0(1,:));
m_contourf(Plg,Plt,data'./4.3.*100000,20,'linestyle','none');shading flat;%colormahold on
shading flat;colormap(m_colmap('divergent'));hold on
hold on;m_contour(Plg,Plt,data',[0 0],'linestyle','--','color','black');
% hold on;m_contour(Plg,Plt,data',[5 10],'linestyle','-','color','black');
% === F test
dataP0 = squeeze(abs(r_matriAllSummerltN2NCEP(:,:)));
% dataP0(land==1) = 0;
dataP = cat(1,dataP0,dataP0(1,:));
[Plg2,Plt2]=meshgrid(lon,latNCEP);
Plt2(dataP'<0.4)=nan;%95% FOR 22
h1=m_plot(Plg2,Plt2,'.','color',[0 .5 0],'Markersize',2);
% m_coast('color','k','linewidth',1.5);
m_coast('patch',[.7 .7 .7],'edgecolor',[.7 .7 .7]);
m_grid('box','fancy','FontSize',8,'FontWeight','bold','xtick',[-180:40:360],...
                'ytick',[30:10:90],'linest',':','tickdir','in');% add a standard colorbar.
caxis([-2 2])
% ===========================================
dataP = cat(1,maskNor,maskNor(1,:));
[x,y] =  find(dataP==6);
kk = boundary(x,y);
% m_plot(lon(x(kk)),lat(y(kk)),'k','linewidth',1)
[x,y] =  find(dataP==1);
kk = boundary(x,y);
m_plot(lon(x(kk)),lat(y(kk)),'y','linewidth',1)
dataP2 = circshift(dataP,70,1);
lon2 = circshift(lon,70,1);
[x,y] =  find(dataP2==3);
kk = boundary(x,y);
m_plot(lon2(x(kk)),lat(y(kk)),'r','linewidth',1)
lon2 = [0:360];
lat2 = lon2.*0+85;
m_plot(lon2,lat2,'Color',[0.6 0.2 1],'linewidth',1)
% ===========================================
% a
text(-0.65,0.65 ,{'a NCEP'},...
    'LineStyle','none',...
    'FontWeight','bold',...
    'FontSize',8,...
    'FontName','Helvetica');
title({'Summer'},'FontSize',8,'position',[0,0.75]);
%
axes(ha(2));
data0 = diffLtN2WinNCEP;
data = cat(1,data0,data0(1,:));
m_contourf(Plg,Plt,data'./4.3.*100000,20,'linestyle','none');shading flat;%colormahold on
shading flat;colormap(m_colmap('divergent'));hold on
hold on;m_contour(Plg,Plt,data',[0 0],'linestyle','--','color','black');
% hold on;m_contour(Plg,Plt,data',[5 10],'linestyle','-','color','black');
% m_coast('color','k','linewidth',1.5);
% m_coast('patch',[.7 .7 .7],'edgecolor',[.7 .7 .7]);
dataP0 = squeeze(abs(r_matriAllWinltN2NCEP(:,:)));
% dataP0(land==1) = 0;
dataP = cat(1,dataP0,dataP0(1,:));
[Plg2,Plt2]=meshgrid(lon,latNCEP);
Plt2(dataP'<0.44)=nan;%95% FOR 22
h1=m_plot(Plg2,Plt2,'.','color',[0 .5 0],'Markersize',2);
m_coast('patch',[.7 .7 .7],'edgecolor',[.7 .7 .7]);
% m_grid('xtick',12,'tickdir','out','ytick',[20 35 50 65 80 90],'fontsize',12,'linewi',1.,'linest','--');
m_grid('box','fancy','FontSize',8,'FontWeight','bold','Fontcolor','white','xtick',[-180:40:360],...
                'ytick',[30:10:90],'linest',':','tickdir','in');% add a standard colorbar.
dataP = cat(1,maskNor,maskNor(1,:));
[x,y] =  find(dataP==6);
kk = boundary(x,y);
% m_plot(lon(x(kk)),lat(y(kk)),'k','linewidth',1)
[x,y] =  find(dataP==1);
kk = boundary(x,y);
m_plot(lon(x(kk)),lat(y(kk)),'y','linewidth',1)
dataP2 = circshift(dataP,70,1);
lon2 = circshift(lon,70,1);
[x,y] =  find(dataP2==3);
kk = boundary(x,y);
m_plot(lon2(x(kk)),lat(y(kk)),'r','linewidth',1)
lon2 = [0:360];
lat2 = lon2.*0+85;
m_plot(lon2,lat2,'Color',[0.6 0.2 1],'linewidth',1)
title({'Winter'},'FontSize',8,'position',[0,0.75],'FontName','Helvetica');
text(-0.65,0.65 ,{'b NCEP'},...
    'LineStyle','none',...
    'FontWeight','bold',...
    'FontSize',8,...
    'FontName','Helvetica');
caxis([-2 2])
%
axes(ha(3));
data0 = diffLtN2SumERA5;
data = cat(1,data0,data0(1,:));
m_contourf(Plg,Plt,data'./4.3.*100000,20,'linestyle','none');shading flat;%colormahold on
shading flat;colormap(m_colmap('divergent'));hold on
hold on;m_contour(Plg,Plt,data',[0 0],'linestyle','--','color','black');
% hold on;m_contour(Plg,Plt,data',[5 10],'linestyle','-','color','black');
% === F test
dataP0 = squeeze(abs(r_matriAllSummerltN2ERA5(:,:)));
% dataP0(land==1) = 0;
dataP = cat(1,dataP0,dataP0(1,:));
[Plg2,Plt2]=meshgrid(lon,latNCEP);
Plt2(dataP'<0.4)=nan;%95% FOR 22
h1=m_plot(Plg2,Plt2,'.','color',[0 .5 0],'Markersize',2);
% m_coast('color','k','linewidth',1.5);
m_coast('patch',[.7 .7 .7],'edgecolor',[.7 .7 .7]);
m_grid('box','fancy','FontSize',8,'FontWeight','bold','xtick',[-180:40:360],...
                'ytick',[30:10:90],'linest',':','tickdir','in');% add a standard colorbar.
caxis([-2 2])
dataP = cat(1,maskNor,maskNor(1,:));
[x,y] =  find(dataP==6);
kk = boundary(x,y);
% m_plot(lon(x(kk)),lat(y(kk)),'k','linewidth',1)
[x,y] =  find(dataP==1);
kk = boundary(x,y);
m_plot(lon(x(kk)),lat(y(kk)),'y','linewidth',1)
dataP2 = circshift(dataP,70,1);
lon2 = circshift(lon,70,1);
[x,y] =  find(dataP2==3);
kk = boundary(x,y);
m_plot(lon2(x(kk)),lat(y(kk)),'r','linewidth',1)
lon2 = [0:360];
lat2 = lon2.*0+85;
m_plot(lon2,lat2,'Color',[0.6 0.2 1],'linewidth',1)
text(-0.65,0.65 ,{'c ERA5'},...
    'LineStyle','none',...
    'FontWeight','bold',...
    'FontSize',8,...
    'FontName','Helvetica');
%
data0 = diffLtN2WinERA5;
data = cat(1,data0,data0(1,:));
axes(ha(4));
m_contourf(Plg,Plt,data'./4.3.*100000,20,'linestyle','none');shading flat;%colormahold on
shading flat;colormap(m_colmap('divergent'));hold on
hold on;m_contour(Plg,Plt,data',[0 0],'linestyle','--','color','black');
% hold on;m_contour(Plg,Plt,data',[5 10],'linestyle','-','color','black');
m_coast('color','k','linewidth',1.5);
dataP0 = squeeze(abs(r_matriAllWinltN2ERA5(:,:)));
% dataP0(land==1) = 0;
dataP = cat(1,dataP0,dataP0(1,:));
[Plg2,Plt2]=meshgrid(lon,latNCEP);
Plt2(dataP'<0.44)=nan;%95% FOR 22
h1=m_plot(Plg2,Plt2,'.','color',[0 .5 0],'Markersize',2);
% m_coast('color','k','linewidth',1.5);
m_coast('patch',[.7 .7 .7],'edgecolor',[.7 .7 .7]);
% m_grid('xtick',12,'tickdir','out','ytick',[20 35 50 65 80 90],'fontsize',12,'linewi',1.,'linest','--');
m_grid('box','fancy','FontSize',8,'FontWeight','bold','Fontcolor','white','xtick',[-180:40:360],...
                'ytick',[30:10:90],'linest',':','tickdir','in');% add a standard colorbar.
dataP = cat(1,maskNor,maskNor(1,:));
[x,y] =  find(dataP==6);
kk = boundary(x,y);
% m_plot(lon(x(kk)),lat(y(kk)),'k','linewidth',1)
[x,y] =  find(dataP==1);
kk = boundary(x,y);
m_plot(lon(x(kk)),lat(y(kk)),'y','linewidth',1)
dataP2 = circshift(dataP,70,1);
lon2 = circshift(lon,70,1);
[x,y] =  find(dataP2==3);
kk = boundary(x,y);
m_plot(lon2(x(kk)),lat(y(kk)),'r','linewidth',1)
lon2 = [0:360];
lat2 = lon2.*0+85;
m_plot(lon2,lat2,'Color',[0.6 0.2 1],'linewidth',1)
text(-0.65,0.65 ,{'d ERA5'},...
    'LineStyle','none',...
    'FontWeight','bold',...
    'FontSize',8,...
    'FontName','Helvetica');
caxis([-2 2])

axes(ha(5));
data0 = diffLtN2SumJRA;
data = cat(1,data0,data0(1,:));
m_contourf(Plg,Plt,data'./4.3.*100000,20,'linestyle','none');shading flat;%colormahold on
shading flat;colormap(m_colmap('divergent'));hold on
hold on;m_contour(Plg,Plt,data',[0 0],'linestyle','--','color','black');
% hold on;m_contour(Plg,Plt,data',[5 10],'linestyle','-','color','black');
% === F test
dataP0 = squeeze(abs(r_matriAllSummerltN2JRA(:,:)));
% dataP0(land==1) = 0;
dataP = cat(1,dataP0,dataP0(1,:));
[Plg2,Plt2]=meshgrid(lon,latNCEP);
Plt2(dataP'<0.4)=nan;%95% FOR 22
h1=m_plot(Plg2,Plt2,'.','color',[0 .5 0],'Markersize',2);
% m_coast('color','k','linewidth',1.5);
m_coast('patch',[.7 .7 .7],'edgecolor',[.7 .7 .7]);
m_grid('box','fancy','FontSize',8,'FontWeight','bold','xtick',[-180:40:360],...
                'ytick',[30:10:90],'linest',':','tickdir','in');% add a standard colorbar.
caxis([-2 2])
dataP = cat(1,maskNor,maskNor(1,:));
[x,y] =  find(dataP==6);
kk = boundary(x,y);
% m_plot(lon(x(kk)),lat(y(kk)),'k','linewidth',1)
[x,y] =  find(dataP==1);
kk = boundary(x,y);
m_plot(lon(x(kk)),lat(y(kk)),'y','linewidth',1)
dataP2 = circshift(dataP,70,1);
lon2 = circshift(lon,70,1);
[x,y] =  find(dataP2==3);
kk = boundary(x,y);
m_plot(lon2(x(kk)),lat(y(kk)),'r','linewidth',1)
lon2 = [0:360];
lat2 = lon2.*0+85;
m_plot(lon2,lat2,'Color',[0.6 0.2 1],'linewidth',1)
text(-0.65,0.65 ,{'e JRA55'},...
    'LineStyle','none',...
    'FontWeight','bold',...
    'FontSize',8,...
    'FontName','Helvetica');
%
data0 = diffLtN2WinJRA;
data = cat(1,data0,data0(1,:));
axes(ha(6));
m_contourf(Plg,Plt,data'./4.3.*100000,20,'linestyle','none');shading flat;%colormahold on
shading flat;colormap(m_colmap('divergent'));hold on
hold on;m_contour(Plg,Plt,data',[0 0],'linestyle','--','color','black');
% hold on;m_contour(Plg,Plt,data',[5 10],'linestyle','-','color','black');
% m_coast('color','k','linewidth',1.5);
dataP0 = squeeze(abs(r_matriAllWinltN2JRA(:,:)));
% dataP0(land==1) = 0;
dataP = cat(1,dataP0,dataP0(1,:));
[Plg2,Plt2]=meshgrid(lon,latNCEP);
Plt2(dataP'<0.44)=nan;%95% FOR 22
h1=m_plot(Plg2,Plt2,'.','color',[0 .5 0],'Markersize',2);
m_coast('patch',[.7 .7 .7],'edgecolor',[.7 .7 .7]);
% m_grid('xtick',12,'tickdir','out','ytick',[20 35 50 65 80 90],'fontsize',12,'linewi',1.,'linest','--');
m_grid('box','fancy','FontSize',8,'FontWeight','bold','Fontcolor','white','xtick',[-180:40:360],...
                'ytick',[30:10:90],'linest',':','tickdir','in');% add a standard colorbar.
dataP = cat(1,maskNor,maskNor(1,:));
[x,y] =  find(dataP==6);
kk = boundary(x,y);
% m_plot(lon(x(kk)),lat(y(kk)),'k','linewidth',1)
[x,y] =  find(dataP==1);
kk = boundary(x,y);
m_plot(lon(x(kk)),lat(y(kk)),'y','linewidth',1)
dataP2 = circshift(dataP,70,1);
lon2 = circshift(lon,70,1);
[x,y] =  find(dataP2==3);
kk = boundary(x,y);
m_plot(lon2(x(kk)),lat(y(kk)),'r','linewidth',1)
lon2 = [0:360];
lat2 = lon2.*0+85;
m_plot(lon2,lat2,'Color',[0.6 0.2 1],'linewidth',1)
text(-0.65,0.65 ,{'f JRA55'},...
    'LineStyle','none',...
    'FontWeight','bold',...
    'FontSize',8,...
    'FontName','Helvetica');
caxis([-2 2])

h2=colorbar('v');%horizontal
set(h2,'position',[0.9 0.08 0.02 0.85]); % ??position:?left??bottom??width??height
set(gca,'FontSize',8,'FontWeight','bold','LineWidth',1.5,'FontName','Helvetica');
set(get(h2,'label'),'string','10^{-5} s^{-2}/decade','FontWeight','bold',...
    'rotation',0,'position',[1.,2.13])
% sgtitle(['Fig. 2 The trend in the static stability of the lower atmosphere in the Arctic during 1979–2021.'],...
%     'fontangle','normal','FontSize',8,'FontName','Helvetica')
sav_name=['npj/final submission/FIG2trendNew'];
print ('-r600','-depsc2',(sav_name))
print ('-r600','-dpdf',(sav_name))
print ('-r600','-djpeg',(sav_name))
%% =====
% ==== pick up the N2 points for region-A, region-B and
% polar
% ==== for mask3 maskNor=3. Norweigh sea
yrs = 42;
lmorAllNorWinNCEP = ones(xmJRA*17,yrs).*nan;
lmorAllNorSumNCEP = ones(xmJRA*17,yrs).*nan;
lmorAllIceWinNCEP = ones(xmJRA*17,yrs).*nan;
lmorAllIceSumNCEP = ones(xmJRA*17,yrs).*nan;
lmorAllPolWinNCEP = ones(xmJRA*17,yrs).*nan;
lmorAllPolSumNCEP = ones(xmJRA*17,yrs).*nan;
lmorAllBriWinNCEP = ones(xmJRA*17,yrs).*nan;
lmorAllBriSumNCEP = ones(xmJRA*17,yrs).*nan;

lmorAllNorWinERA5 = ones(xmJRA*17,yrs).*nan;
lmorAllNorSumERA5 = ones(xmJRA*17,yrs).*nan;
lmorAllIceWinERA5 = ones(xmJRA*17,yrs).*nan;
lmorAllIceSumERA5 = ones(xmJRA*17,yrs).*nan;
lmorAllPolWinERA5 = ones(xmJRA*17,yrs).*nan;
lmorAllPolSumERA5 = ones(xmJRA*17,yrs).*nan;
lmorAllBriWinERA5 = ones(xmJRA*17,yrs).*nan;
lmorAllBriSumERA5 = ones(xmJRA*17,yrs).*nan;

lmorAllNorWinJRA = ones(xmJRA*17,yrs).*nan;
lmorAllNorSumJRA = ones(xmJRA*17,yrs).*nan;
lmorAllIceWinJRA = ones(xmJRA*17,yrs).*nan;
lmorAllIceSumJRA = ones(xmJRA*17,yrs).*nan;
lmorAllPolWinJRA = ones(xmJRA*17,yrs).*nan;
lmorAllPolSumJRA = ones(xmJRA*17,yrs).*nan;
lmorAllBriWinJRA = ones(xmJRA*17,yrs).*nan;
lmorAllBriSumJRA = ones(xmJRA*17,yrs).*nan;

kk1 = 1;
kk2 = 1;
kk3 = 1;
for i = 1:xmJRA
    for j = 1:25
        % === calculate the average amplitude without the trend
        ln2mWinNCEP = squeeze(WinLtN20NCEP(i,j,1:yrs));
        ln2mSumNCEP = squeeze(SummerLtN20NCEP(i,j,1:yrs));
        ln2mWinERA5 = squeeze(WinLtN20ERA5(i,j,1:yrs));
        ln2mSumERA5 = squeeze(SummerLtN20ERA5(i,j,1:yrs));
        ln2mWinJRA = squeeze(WinLtN20JRA(i,j,1:yrs));
        ln2mSumJRA = squeeze(SummerLtN20JRA(i,j,1:yrs));
        if (maskNor(i,j)==3)
            lmorAllNorWinNCEP(kk1,:) = ln2mWinNCEP;
            lmorAllNorSumNCEP(kk1,:) = ln2mSumNCEP;
            lmorAllNorWinERA5(kk1,:) = ln2mWinERA5;
            lmorAllNorSumERA5(kk1,:) = ln2mSumERA5;
            lmorAllNorWinJRA(kk1,:) = ln2mWinJRA;
            lmorAllNorSumJRA(kk1,:) = ln2mSumJRA;
            kk1 = kk1+1;
        elseif (maskNor(i,j)==1)
            lmorAllIceWinNCEP(kk2,:) = ln2mWinNCEP;
            lmorAllIceSumNCEP(kk2,:) = ln2mSumNCEP;
            lmorAllIceWinERA5(kk2,:) = ln2mWinERA5;
            lmorAllIceSumERA5(kk2,:) = ln2mSumERA5;
            lmorAllIceWinJRA(kk2,:) = ln2mWinJRA;
            lmorAllIceSumJRA(kk2,:) = ln2mSumJRA;
            kk2 = kk2+1;        
        elseif (maskNor(i,j)==4)
            lmorAllPolWinNCEP(kk3,:) = ln2mWinNCEP;
            lmorAllPolSumNCEP(kk3,:) = ln2mSumNCEP;
            lmorAllPolWinERA5(kk3,:) = ln2mWinERA5;
            lmorAllPolSumERA5(kk3,:) = ln2mSumERA5;
            lmorAllPolWinJRA(kk3,:) = ln2mWinJRA;
            lmorAllPolSumJRA(kk3,:) = ln2mSumJRA;
            kk3 = kk3+1;
        end
    end
end
%%
% === fig 4 ocsillation N2
figure;
Nh=3; Nw=2; gap=0.05; marg_h=0.08; marg_w=0.05;
% ha = tight_subplot(Nh, Nw, gap, marg_h, marg_w);%
ha = tight_subplot(Nh, Nw, gap, [0.06,0.08], [0.08,0.02]);%
% === the Summer 
x =[1:42];
axes(ha(1));
ydata = nanmean(lmorAllPolSumNCEP./nanmean(lmorAllPolSumNCEP,'all')).*100;
aa= polyfit(x,ydata,1);
plot(x,x.*0,'k','linewidth',1);hold on
test = detrend(lmorAllPolSumNCEP')';
lmorNCEP = test;%./nanmean(lmorAllPolSumNCEP).*100;%
lmmeanPolNCEP = nanmean(lmorNCEP);%regional average
hold on; plot(lmmeanPolNCEP,'r','linewidth',1.5)
[pks,~]=findpeaks(lmmeanPolNCEP);
mean(pks)
[pks,~]=findpeaks(-lmmeanPolNCEP);
mean(pks)
%
test = detrend(lmorAllPolSumERA5')';
lmorERA5 = test;%./nanmean(lmorAllPolSumERA5).*100;%
lmmeanPolERA5 = nanmean(lmorERA5);%regional average
hold on; plot(lmmeanPolERA5,'b','linewidth',1.5)
%
test = detrend(lmorAllPolSumJRA')';
lmorJRA = test;%./nanmean(lmorAllPolSumJRA).*100;%
lmmeanPolJRA = nanmean(lmorJRA);%regional average
hold on; plot(lmmeanPolJRA,'k','linewidth',1.5)
%
set(gca,'FontName','Helvetica','FontSize',8,'LineWidth',1.,'TickDir','in',...
    'XMinorTick','on','TickLength',[0.025 0.05],'XTick',[2:5:42],'XTickLabel',...
    {'1980','1985','1990','1995','2000','2005','2010','2015','2020'},...
    'YMinorTick','on');
xlim([1 42])
ylabel({'N^2 oscillation(s^{-2})'})
gtext({'Polar Region'},'FontSize',8,'FontName','Helvetica');

gtext({'a'},'FontSize',8,'FontWeight','bold','FontName','Helvetica');

axes(ha(3));
ydata = nanmean(lmorAllNorSumNCEP./nanmean(lmorAllNorSumNCEP,'all')).*100;
aa= polyfit(x,ydata,1);
plot(x,x.*0,'k','linewidth',1);hold on
test = detrend(lmorAllNorSumNCEP')';
lmorNCEP = test;%./nanmean(lmorAllNorSumNCEP).*100;%
lmmeanNorNCEP = nanmean(lmorNCEP);%regional average
hold on; plot(lmmeanNorNCEP,'r','linewidth',1.5)
%
test = detrend(lmorAllNorSumERA5')';
lmorERA5 = test;%./nanmean(lmorAllNorSumERA5).*100;%
lmmeanNorERA5 = nanmean(lmorERA5);%regional average
hold on; plot(lmmeanNorERA5,'b','linewidth',1.5)
%
test = detrend(lmorAllNorSumJRA')';
lmorJRA = test;%./nanmean(lmorAllNorSumJRA).*100;%
lmmeanNorJRA = nanmean(lmorJRA);%regional average
hold on; plot(lmmeanNorJRA,'k','linewidth',1.5)
%
set(gca,'FontName','Helvetica','FontSize',8,'LineWidth',1.,'TickDir','in',...
    'XMinorTick','on','TickLength',[0.025 0.05],'XTick',[2:5:42],'XTickLabel',...
    {'1980','1985','1990','1995','2000','2005','2010','2015','2020'},...
    'YMinorTick','on');
xlim([1 42])
ylabel({'N^2 oscillation(s^{-2})'})
gtext({'Region-A'},'FontSize',8,'FontName','Helvetica');

gtext({'c'},'FontSize',8,'FontWeight','bold','FontName','Helvetica');

axes(ha(5));
ydata = nanmean(lmorAllIceSumNCEP./nanmean(lmorAllIceSumNCEP,'all')).*100;
aa= polyfit(x,ydata,1);
plot(x,x.*0,'k','linewidth',1);hold on
test = detrend(lmorAllIceSumNCEP')';
lmorNCEP = test;%./nanmean(lmorAllIceSumNCEP).*100;%
lmmeanIceNCEP = nanmean(lmorNCEP);%regional average
hold on; plot(lmmeanIceNCEP,'r','linewidth',1.5)
%
test = detrend(lmorAllIceSumERA5')';
lmorERA5 = test;%./nanmean(lmorAllIceSumERA5).*100;%
lmmeanIceERA5 = nanmean(lmorERA5);%regional average
hold on; plot(lmmeanIceERA5,'b','linewidth',1.5)
%
test = detrend(lmorAllIceSumJRA')';
lmorJRA = test;%./nanmean(lmorAllIceSumJRA).*100;%
lmmeanIceJRA = nanmean(lmorJRA);%regional average
hold on; plot(lmmeanIceJRA,'k','linewidth',1.5)
set(gca,'FontName','Helvetica','FontSize',8,'LineWidth',1.,'TickDir','in',...
    'XMinorTick','on','TickLength',[0.025 0.05],'XTick',[2:5:42],'XTickLabel',...
    {'1980','1985','1990','1995','2000','2005','2010','2015','2020'},...
    'YMinorTick','on');
xlim([1 42])
legend('0-line','NCEP','ERA5','JRA55','Location','best','NumColumns',2)
ylabel({'N^2 oscillation(s^{-2})'})
xlabel('Time(Year)')
gtext({'Region-B'},'FontSize',8,'FontName','Helvetica');

gtext({'e'},'FontSize',8,'FontWeight','bold','FontName','Helvetica');

% === the winter 
axes(ha(2));
ydata = nanmean(lmorAllPolWinNCEP./nanmean(lmorAllPolWinNCEP,'all')).*100;
aa= polyfit(x,ydata,1);
plot(x,x.*0,'k','linewidth',1);hold on
test = detrend(lmorAllPolWinNCEP')';
lmorNCEP = test;%./nanmean(lmorAllPolWinNCEP).*100;%
lmmeanPolNCEP = nanmean(lmorNCEP);%regional average
hold on; plot(lmmeanPolNCEP,'r','linewidth',1.5)
[pks,~]=findpeaks(lmmeanPolNCEP);
mean(pks)
std(pks)
[pks,~]=findpeaks(-lmmeanPolNCEP);
mean(pks)
std(pks)
%
test = detrend(lmorAllPolWinERA5')';
lmorERA5 = test;%./nanmean(lmorAllPolWinERA5).*100;%
lmmeanPolERA5 = nanmean(lmorERA5);%regional average
hold on; plot(lmmeanPolERA5,'b','linewidth',1.5)
%
test = detrend(lmorAllPolWinJRA')';
lmorJRA = test;%./nanmean(lmorAllPolWinJRA).*100;%
lmmeanPolJRA = nanmean(lmorJRA);%regional average
hold on; plot(lmmeanPolJRA,'k','linewidth',1.5)
%
set(gca,'FontName','Helvetica','FontSize',8,'LineWidth',1.,'TickDir','in',...
    'XMinorTick','on','TickLength',[0.025 0.05],'XTick',[1:5:42],'XTickLabel',...
    {'1980','1985','1990','1995','2000','2005','2010','2015','2020'},...
    'YMinorTick','on');
xlim([1 42])
% ylabel({'N^2 oscillation(s^{-2})'})
% gtext({['Winter']},'FontName','Helvetica','FontSize',7
gtext({'Polar Region'},'FontSize',8,'FontName','Helvetica');

gtext({'b'},'FontSize',8,'FontWeight','bold','FontName','Helvetica');

axes(ha(4));
ydata = nanmean(lmorAllNorWinNCEP./nanmean(lmorAllNorWinNCEP,'all')).*100;
aa= polyfit(x,ydata,1);
plot(x,x.*0,'k','linewidth',1);hold on
test = detrend(lmorAllNorWinNCEP')';
lmorNCEP = test;%./nanmean(lmorAllNorWinNCEP).*100;%
lmmeanNorNCEP = nanmean(lmorNCEP);%regional average
hold on; plot(lmmeanNorNCEP,'r','linewidth',1.5)
%
test = detrend(lmorAllNorWinERA5')';
lmorERA5 = test;%./nanmean(lmorAllNorWinERA5).*100;%
lmmeanNorERA5 = nanmean(lmorERA5);%regional average
hold on; plot(lmmeanNorERA5,'b','linewidth',1.5)
%
test = detrend(lmorAllNorWinJRA')';
lmorJRA = test;%./nanmean(lmorAllNorWinJRA).*100;%
lmmeanNorJRA = nanmean(lmorJRA);%regional average
hold on; plot(lmmeanNorJRA,'k','linewidth',1.5)
%
set(gca,'FontName','Helvetica','FontSize',8,'LineWidth',1.,'TickDir','in',...
    'XMinorTick','on','TickLength',[0.025 0.05],'XTick',[1:5:42],'XTickLabel',...
    {'1980','1985','1990','1995','2000','2005','2010','2015','2020'},...
    'YMinorTick','on');
xlim([1 42])
% ylabel({'N^2 oscillation(s^{-2})'})
gtext({'Region-A'},'FontSize',8,'FontName','Helvetica');

gtext({'d'},'FontSize',8,'FontWeight','bold','FontName','Helvetica');

axes(ha(6));
ydata = nanmean(lmorAllIceWinNCEP./nanmean(lmorAllIceWinNCEP,'all')).*100;
aa= polyfit(x,ydata,1);
plot(x,x.*0,'k','linewidth',1);hold on
test = detrend(lmorAllIceWinNCEP')';
lmorNCEP = test;%./nanmean(lmorAllIceWinNCEP).*100;%
lmmeanIceNCEP = nanmean(lmorNCEP);%regional average
hold on; plot(lmmeanIceNCEP,'r','linewidth',1.5)
%
test = detrend(lmorAllIceWinERA5')';
lmorERA5 = test;%./nanmean(lmorAllIceWinERA5).*100;%
lmmeanIceERA5 = nanmean(lmorERA5);%regional average
hold on; plot(lmmeanIceERA5,'b','linewidth',1.5)
%
test = detrend(lmorAllIceWinJRA')';
lmorJRA = test;%./nanmean(lmorAllIceWinJRA).*100;%
lmmeanIceJRA = nanmean(lmorJRA);%regional average
hold on; plot(lmmeanIceJRA,'k','linewidth',1.5)
set(gca,'FontName','Helvetica','FontSize',8,'LineWidth',1.,'TickDir','in',...
    'XMinorTick','on','TickLength',[0.025 0.05],'XTick',[1:5:42],'XTickLabel',...
    {'1980','1985','1990','1995','2000','2005','2010','2015','2020'},...
    'YMinorTick','on');
xlim([1 42])
% ylabel({'N^2 oscillation(s^{-2})'})
xlabel('Time(Year)')
gtext({'Region-B'},'FontSize',8,'FontName','Helvetica');

gtext({'f'},'FontSize',8,'FontWeight','bold','FontName','Helvetica');

set(gcf,'color','white','Units','centimeters','Position',[3 3 18 12]);
% sav_name=['DAY47/LinePlot3RegionRelativeN2_LowTroOsciWinAndSum'];
% sgtitle(['Fig. 4 The relative oscillatory variations in the static stability of the lower atmosphere over three regions.'],...
%     'fontangle','normal','FontSize',8,'FontName','Helvetica')

sav_name=['npj/final submission/FIG4new'];%

print ('-r600','-djpeg',(sav_name))
print ('-r600','-depsc2',(sav_name))
print ('-r600','-dpdf',(sav_name))
%% === fig 6 calculate the ccww500 between AOI and OMEGA 500HPA 

fid = 'NCEPdata\omega.mon.mean.nc';
ncdisp(fid)
time = ncread(fid,'time');
level = ncread(fid,'level');
omega =  double(ncread(fid,'omega'));
[xm,ym,zm,tm] = size(omega);
% - average
momega = zeros(xm,ym,zm,12);
domega = zeros(xm,ym,zm,tm);
for jt =1:12
    momega(:,:,:,jt) = mean(omega(:,:,:,jt:12:360),4);
    domega(:,:,:,jt:12:end) = omega(:,:,:,jt:12:end)-momega(:,:,:,jt);
end
% running average
% 500HpA WW AND AOI cc
nd = 6;
nd2 = 13;
tid = 1;
AOIorRA0 = running_average(AOIor,nd,nd,nd2);
AOIorRA = AOIorRA0(nd+1:tmERA5-nd);
% txx = aoi(nd+1:tm-nd,1)+aoi(nd+1:tm-nd,2)./12;
txx = aoi(nd+tid:tmERA5-nd,1)+aoi(nd+tid:tmERA5-nd,2)./12;
% === the ccww500 between aoi and N2
txx = aoi(nd+tid:tmERA5-nd,1)+aoi(nd+tid:tmERA5-nd,2)./12;
tm05h = length(txx);
% === the ccww500 between aoi and N2
ccww500AOIRANCEP = zeros(xmNCEP,ymNCEP);
ccww500AOINCEP = zeros(xmNCEP,ymNCEP);
for i = 1:xmNCEP
    for j = 1:ymNCEP
        %NCEP
        ln2mAll2 = detrend(squeeze(domega(i,j,6,tid:tmERA5)));
        ln2mAll2RA0 = running_average(ln2mAll2,nd,nd,nd2);
        ln2mAll2RA = ln2mAll2RA0(nd+1:tmERA5-nd);
        ccww500AOINCEP(i,j) = corr(ln2mAll2,AOIor);
        ccww500AOIRANCEP(i,j) = corr(ln2mAll2RA',AOIorRA');
    end
end
%
%% ====== FIG6 a====
figure;
set(gcf,'color','white','Units','centimeters','Position',[3 3 6 6.5]);
Nh=1; Nw=1; gap=0.04; 
ha = tight_subplot(Nh, Nw, [0.04,0.02], [0.17,0.05], [0.1,0.1]);%
m_proj('stereographic','lat',90,'long',-10,'radius',35);
lon = cat(1,lon0NCEP,lon0NCEP(1,:));
[Plg,Plt]=meshgrid(lon,latNCEP);
axes(ha(1));%
data0 = ccww500AOIRANCEP;
data = cat(1,data0,data0(1,:));
m_contourf(Plg,Plt,data',20,'linestyle','none');shading flat;%colormahold on
shading flat;colormap(m_colmap('divergent'));hold on
hold on;m_contour(Plg,Plt,data',[0 0],'linestyle','--','color','black');
hold on;m_contour(Plg,Plt,data',[-0.4227 -0.4227],'linestyle','-','color','black');
m_coast('patch',[.7 .7 .7],'edgecolor',[.7 .7 .7]);
m_grid('box','fancy','fontsize',8,'FontWeight','bold','xtick',[-180:40:360],...
                'ytick',[30:10:90],'linest',':','tickdir','in');% add a standard colorbar.
caxis([-1 1])
% ===========================================
dataP = cat(1,maskNor,maskNor(1,:));
[x,y] =  find(dataP==6);
kk = boundary(x,y);
[x,y] =  find(dataP==1);
kk = boundary(x,y);
m_plot(lon(x(kk)),lat(y(kk)),'y','linewidth',1)
[x,y] =  find(dataP==4);
kk = boundary(x,y);
m_plot(lon(x(kk(149:270))),lat(y(kk(149:270))),'Color',[0.6 0.2 1],'linewidth',1)
dataP2 = circshift(dataP,70,1);
lon2 = circshift(lon,70,1);
[x,y] =  find(dataP2==3);
kk = boundary(x,y);
m_plot(lon2(x(kk)),lat(y(kk)),'r','linewidth',1)
% ===========================================
% a
text(-0.6,0.6 ,{'a'},...
    'LineStyle','none',...
    'FontWeight','bold',...
    'FontSize',8,...
    'FontName','Helvetica');
%CC of N^2 and surface air temperature
h2=colorbar('h');%horizontal
set(h2,'position',[0.1 0.07 0.8 0.04]); % ??position:?left??bottom??width??height
% title({'Summer'},'fontsize',8,'position',[0,0.67]);
sav_name=['npj/final submission/Fig6aCCwwAOINCEP'];
print ('-r600','-djpeg',(sav_name))
print ('-r600','-depsc2',(sav_name))
print ('-r600','-dpdf',(sav_name))


%% === fig6 bc
HCRWW500NCEP = ones(xmJRA*7,tmERA5).*nan;
kk3 = 1;
for i = 1:xmJRA
    for j = 1:25
        % === calculate the average amplitude without the trend
        ln1mNCEP = squeeze(domega(i,j,6,1:tmERA5));
        if (ccww500AOIRANCEP(i,j)<=-0.4227)
            HCRWW500NCEP(kk3,:) = ln1mNCEP;
             kk3 = kk3+1;
        end
    end
end
% ===== fig 6 bc AOI and HCR region WW 500hpa rcc
figure;
set(gcf,'color','white','Units','centimeters','Position',[3 3 6 6.5]);
Nh=2; Nw=1; gap=0.04; 
ha = tight_subplot(Nh, Nw, [0.02,0.01], [0.178,0.06], [0.189,0.13]);%
axes(ha(1));
lm1mean = nanmean(HCRWW500NCEP);
lm2mean = AOIor(1:tmERA5);
lm1mid = running_average(lm1mean,nd,nd,nd2); 
lm2mid = running_average(lm2mean,nd,nd,nd2); 
lm1 = lm1mid(nd+1:tmERA5-nd);
lm2 = lm2mid(nd+1:tmERA5-nd);
[hAx,hLine1,hLine2] = plotyy(txx,lm1,txx,lm2);
hLine1.LineWidth = 1.5;
hLine2.LineWidth = 1.5;
xlim(hAx(1),[1980,2021])
xlim(hAx(2),[1980,2021])
% ylim(hAx(1),[-0.00002 0.00002])
% ylim(hAx(2),[-0.02,0.02])
ylabel(hAx(1),{'omega anomaly(Pa/s)'})
ylabel(hAx(2),'AOI')
set(hAx(1),'FontSize',8,'LineWidth',1,'TickDir',...
    'in','XTick',[1980:5:2021],'XMinorTick','on','XTickLabel',...
    [],'TickLength',[0.025 0.05],...
      'YMinorTick','on','FontName','Helvetica');
set(hAx(2),'FontSize',8,'LineWidth',1,'TickDir',...
    'in','XTick',[1980:5:2021],'XMinorTick','on','TickLength',[0.025 0.05],'XTickLabel',...
    [],...
    'YMinorTick','on');
text('Parent',gca,'VerticalAlignment','baseline',...
    'FontSize',8,'FontName','Helvetica',...
    'String',{'Polar Region'},...
    'Position',[1980.7207729383 0.013]);
text('Parent',gca,'VerticalAlignment','baseline','FontWeight','bold',...
    'FontSize',8,'FontName','Helvetica',...
    'String',{'b'},...
    'Position',[1981.20891518602 -1.8*0.01 0]);
hold on;plot(txx,txx*0,'k-')
%lm1,txx,lm2
axes(ha(2));
nd = 6;
tmm = length(lm2);
syCC2 = zeros(tmm,1);
a1 = mean(lm1(nd:tmm-nd));
a2 = mean(lm2(nd:tmm-nd));
for kk = nd+1:tmm-nd-1
    xk = lm1(kk-nd:kk+nd);
    yk = lm2(kk-nd:kk+nd);
%     a1 = mean(xk);
%     a2 = mean(yk);
    dxk = xk-a1;
    dyk = yk-a2;
    ntr = sum(dxk.*dyk);
    dtr = (sum(dxk.*dxk).*sum(dyk.*dyk))^0.5;
    syCC2(kk) = ntr./dtr;
end
area(txx,syCC2)
hold on;plot(txx,txx*0+0.35,'k--')
hold on;plot(txx,txx*0-0.35,'k--')
set(gca,'FontSize',8,'LineWidth',1,'FontName','Helvetica');
set(gca,'FontName','Helvetica','FontSize',8,'TickDir','in',...
    'XMinorTick','on','YMinorTick','on','TickLength',[0.025 0.05],'XTick',[1980:5:2021],'XTickLabel',...
    {'1980','1985','1990','1995','2000','2005','2010','2015','2020'});
ylabel({'Correlation ','coefficient'})
xlabel('Year')
xlim([1980,2021])
text('Parent',gca,'VerticalAlignment','baseline','FontWeight','bold',...
    'FontSize',8,'FontName','Helvetica',...
    'String',{'c'},...
    'Position',[1980.83601087494 -0.867539046010974 0]);
sav_name=['npj/final submission/Fig6bcRCCwwAOINCEP'];

print ('-r600','-djpeg',(sav_name))
print ('-r600','-depsc2',(sav_name))
print ('-r600','-dpdf',(sav_name))
%% fig 8
% === pick out the each data in region-a and region-b based on masknor
tm2 = 504-12;
dLTN20NCEP_AA = ones(xmNCEP*7,tmERA5).*nan;
dLTN20NCEP_BB = ones(xmNCEP*7,tmERA5).*nan;
kk1 = 1;
kk2 = 1;
for i = 1:xmNCEP
    for j = 1:25
        % === calculate the average amplitude without the trend
        ln1mNCEP = squeeze(dLTN20NCEP(i,j,1:tmERA5));
        if (maskNor(i,j)==3)%norwegian
            dLTN20NCEP_AA(kk1,:) = ln1mNCEP;
            kk1 = kk1+1;
        elseif (maskNor(i,j)==1)%iceland
            dLTN20NCEP_BB(kk2,:) = ln1mNCEP;
            kk2 = kk2+1;
        end
    end
end
% cc distribution of N2 in each point and region-a  or region-b

lm1Ag = nanmean(dLTN20NCEP_AA);
lm2Ag = nanmean(dLTN20NCEP_BB);
lm1mid = running_average(lm1Ag,nd,nd,nd2); 
lm2mid = running_average(lm2Ag,nd,nd,nd2); 
lm1 = lm1mid(nd+1:tmERA5-nd);
lm2 = lm2mid(nd+1:tmERA5-nd);
CCn2AA = zeros(xmNCEP,ymNCEP);
CCn2BB = zeros(xmNCEP,ymNCEP);
for i = 1:xmNCEP
    for j = 1:25
        lm3Ag = squeeze(dLTN20NCEP(i,j,:));
        lm3mid = running_average(lm3Ag,nd,nd,nd2); 
        lm3 = lm3mid(nd+1:tmERA5-nd);
        CCn2AA(i,j) = corr(lm1',lm3');
        CCn2BB(i,j) = corr(lm2',lm3');       
    end
end
%%
figure;
set(gcf,'color','white','Units','centimeters','Position',[3 3 12 6]);
Nh=1; Nw=2; gap=0.04; 
ha = tight_subplot(Nh, Nw, [0.06,0.02], [0.08,0.08], [0.025,0.1]);%
m_proj('stereographic','lat',90,'long',-10,'radius',60);
lon = cat(1,lon0NCEP,lon0NCEP(1,:));
[Plg,Plt]=meshgrid(lon,latNCEP);
axes(ha(1));%
data0 = CCn2AA;
data = cat(1,data0,data0(1,:));
m_contourf(Plg,Plt,data',20,'linestyle','none');shading flat;%colormahold on
shading flat;colormap(m_colmap('divergent'));hold on
hold on;m_contour(Plg,Plt,data',[0 0],'linestyle','--','color','black');
hold on;m_contour(Plg,Plt,data',[-0.4227 0.4227],'linestyle','-','color','black');
% m_coast('color','k','linewidth',1.5);
m_coast('patch',[.7 .7 .7],'edgecolor',[.7 .7 .7]);
m_grid('box','fancy','FontSize',8,'FontWeight','bold','xtick',[-180:40:360],...
                'ytick',[30:10:90],'linest',':','tickdir','in');% add a standard colorbar.
caxis([-1 1])
% ===========================================
dataP = cat(1,maskNor,maskNor(1,:));
[x,y] =  find(dataP==6);
kk = boundary(x,y);
[x,y] =  find(dataP==1);
kk = boundary(x,y);
m_plot(lon(x(kk)),lat(y(kk)),'y','linewidth',1)
[x,y] =  find(dataP==4);
kk = boundary(x,y);
m_plot(lon(x(kk(149:270))),lat(y(kk(149:270))),'Color',[0.6 0.2 1],'linewidth',1)
dataP2 = circshift(dataP,70,1);
lon2 = circshift(lon,70,1);
[x,y] =  find(dataP2==3);
kk = boundary(x,y);
m_plot(lon2(x(kk)),lat(y(kk)),'r','linewidth',1)
% ===========================================
% a
text(-1.15,1.15 ,{'a'},...
    'LineStyle','none',...
    'FontWeight','bold',...
    'FontSize',8,...
    'FontName','Helvetica');
axes(ha(2));%
data0 = CCn2BB;
data = cat(1,data0,data0(1,:));
m_contourf(Plg,Plt,data',20,'linestyle','none');shading flat;%colormahold on
shading flat;colormap(m_colmap('divergent'));hold on
hold on;m_contour(Plg,Plt,data',[0 0],'linestyle','--','color','black');
hold on;m_contour(Plg,Plt,data',[-0.4227 0.4227],'linestyle','-','color','black');
% m_coast('color','k','linewidth',1.5);
m_coast('patch',[.7 .7 .7],'edgecolor',[.7 .7 .7]);
m_grid('box','fancy','FontSize',8,'FontWeight','bold','xtick',[-180:40:360],...
                'ytick',[30:10:90],'linest',':','tickdir','in');% add a standard colorbar.
caxis([-1 1])
% ===========================================
dataP = cat(1,maskNor,maskNor(1,:));
[x,y] =  find(dataP==6);
kk = boundary(x,y);
[x,y] =  find(dataP==1);
kk = boundary(x,y);
m_plot(lon(x(kk)),lat(y(kk)),'y','linewidth',1)
[x,y] =  find(dataP==4);
kk = boundary(x,y);
m_plot(lon(x(kk(149:270))),lat(y(kk(149:270))),'Color',[0.6 0.2 1],'linewidth',1)
dataP2 = circshift(dataP,70,1);
lon2 = circshift(lon,70,1);
[x,y] =  find(dataP2==3);
kk = boundary(x,y);
m_plot(lon2(x(kk)),lat(y(kk)),'r','linewidth',1)
% ===========================================
% a
text(-1.15,1.15 ,{'b'},...
    'LineStyle','none',...
    'FontWeight','bold',...
    'FontSize',8,...
    'FontName','Helvetica');
h2=colorbar;%horizontal
set(h2,'position',[0.923 0.1 0.02 0.8]); % ??position:?left??bottom??width??height
% title({'Summer'},'FontSize',8,'position',[0,0.67]);
sav_name=['npj/final submission/Fig8CCN2NCEP'];
print ('-r600','-djpeg',(sav_name))
print ('-r600','-depsc2',(sav_name))
print ('-r600','-dpdf',(sav_name))
%% === fig 7
% SST\SH\LH plot the 7 year
tm5 = 1826;
LTN20NCEP_AA90_94 = ones(32,tm5).*nan;
LTN20NCEP_BB90_94 = ones(15,tm5).*nan;
SST_AA90_94 = ones(32, tm5).*nan;
SST_BB90_94 = ones(15, tm5).*nan;
air925_AA90_94 = ones(32, tm5).*nan;
air925_BB90_94 = ones(15, tm5).*nan;
air850_AA90_94 = ones(32, tm5).*nan;
air850_BB90_94 = ones(15, tm5).*nan;
omega_AA90_94 = ones(32,tm5).*nan;
omega_BB90_94 = ones(15,tm5).*nan;
SH_AA90_94 = ones(32, tm5).*nan;
SH_BB90_94 = ones(15, tm5).*nan;
LH_AA90_94 = ones(32, tm5).*nan;
LH_BB90_94 = ones(15, tm5).*nan;
UL_AA90_94 = ones(32, tm5).*nan;
DL_AA90_94 = ones(32, tm5).*nan;
US_AA90_94 = ones(32, tm5).*nan;
DS_AA90_94 = ones(32, tm5).*nan;
UL_BB90_94 = ones(15, tm5).*nan;
DL_BB90_94 = ones(15, tm5).*nan;
US_BB90_94 = ones(15, tm5).*nan;
DS_BB90_94 = ones(15, tm5).*nan;
days = [365,365,366,365,365];
    k = 1;
for yr = 1990:1994
    tmm = days(k);
    fnm = ['DATA/resultAABB',num2str(yr),'.mat'];
    load(fnm)
    tid1 = sum(days(1:k))-days(k)+1;
    tid2 = sum(days(1:k));
    LTN20NCEP_AA90_94(:,tid1:tid2) =    LTN20NCEP_AAYR(1:32,1:tmm);
    LTN20NCEP_BB90_94(:,tid1:tid2) =    LTN20NCEP_BBYR(1:15,1:tmm);
    SST_AA90_94(:,tid1:tid2) =     SST_AAYR(1:32,1:tmm);
    SST_BB90_94(:,tid1:tid2) =     SST_BBYR(1:15,1:tmm);
    air925_AA90_94(:,tid1:tid2) =     air925_AAYR(1:32,1:tmm);
    air925_BB90_94(:,tid1:tid2) =     air925_BBYR(1:15,1:tmm);
    air850_AA90_94(:,tid1:tid2) =     air850_AAYR(1:32,1:tmm);
    air850_BB90_94(:,tid1:tid2) =     air850_BBYR(1:15,1:tmm);
    omega_AA90_94(:,tid1:tid2) =    omega_AAYR(1:32,1:tmm);
    omega_BB90_94(:,tid1:tid2) =    omega_BBYR(1:15,1:tmm);
    SH_AA90_94(:,tid1:tid2) =     SH_AAYR(1:32,1:tmm);
    SH_BB90_94(:,tid1:tid2) =     SH_BBYR(1:15,1:tmm);
    LH_AA90_94(:,tid1:tid2) =     LH_AAYR(1:32,1:tmm);
    LH_BB90_94(:,tid1:tid2) =     LH_BBYR(1:15,1:tmm);
    UL_AA90_94(:,tid1:tid2) =     UL_AAYR(1:32,1:tmm);
    DL_AA90_94(:,tid1:tid2) =     DL_AAYR(1:32,1:tmm);
    US_AA90_94(:,tid1:tid2) =     US_AAYR(1:32,1:tmm);
    DS_AA90_94(:,tid1:tid2) =     DS_AAYR(1:32,1:tmm);
    UL_BB90_94(:,tid1:tid2) =     UL_BBYR(1:15,1:tmm);
    DL_BB90_94(:,tid1:tid2) =     DL_BBYR(1:15,1:tmm);
    US_BB90_94(:,tid1:tid2) =     US_BBYR(1:15,1:tmm);
    DS_BB90_94(:,tid1:tid2) =     DS_BBYR(1:15,1:tmm);
    k = k+1;
end
%%
figure;
Nh=3; Nw=2; gap=0.05; marg_h=0.08; marg_w=0.05;
set(gcf,'color','white','Units','centimeters','Position',[3 3 18 12]);
ha = tight_subplot(Nh, Nw, gap, [0.08,0.08], [0.08,0.02]);%
% sst\ul and at lag and head correlation running-average 7-days
axes(ha(1));
lm1Ag = mean(LTN20NCEP_AA90_94,'omitnan');
lm2Ag = mean(SH_AA90_94,'omitnan');
nd = 3;nd2 =7;
lm1mid = running_average(lm1Ag,nd,nd,nd2);
lm2mid = running_average(lm2Ag,nd,nd,nd2);
tmYR = length(lm1mid);
% lm1 = detrend(lm1mid(nd+1:tmYR-nd));
% lm2 = detrend(lm2mid(nd+1:tmYR-nd));
lm1 = (lm1mid(nd+1:tmYR-nd));
lm2 = (lm2mid(nd+1:tmYR-nd));
corrcoef(lm1Ag',lm2Ag')
corrcoef(lm1',lm2')
tmmid = length(lm1);
nd = 14;
l1 = lm1(nd+1:tmmid-2*nd);
lagCC = zeros(2*nd+1,1);
for lg = -nd:nd
    l2 = lm2(nd+1+lg:tmmid-2*nd+lg);
    mid = corrcoef(l1',l2');
    lagCC(lg+nd+1) = mid(1,2);
end
txx2 = [1979.5:1./12:2021.49];
plot([-nd:nd],lagCC,'linewidth',1.5);hold on

ylabel('CC')
set(gca,'FontSize',8,'LineWidth',1,'FontName','Helvetica');
set(gca,'FontName','Helvetica','FontSize',8,'LineWidth',1.,'TickDir','in',...
    'XMinorTick','on','TickLength',[0.025 0.05],...
    'YMinorTick','on');
xlim([-nd,nd])
grid on
gtext('a Region-A SH and N^2','FontSize',8,'FontName','Helvetica','fontweight','bold')

axes(ha(2));
lm1Ag = mean(LTN20NCEP_BB90_94,'omitnan');
lm2Ag = mean(SH_BB90_94,'omitnan');
nd = 3;nd2 =7;
lm1mid = running_average(lm1Ag,nd,nd,nd2);
lm2mid = running_average(lm2Ag,nd,nd,nd2);
% lm1 = detrend(lm1mid(nd+1:tmYR-nd));
% lm2 = detrend(lm2mid(nd+1:tmYR-nd));
lm1 = (lm1mid(nd+1:tmYR-nd));
lm2 = (lm2mid(nd+1:tmYR-nd));
corrcoef(lm1Ag',lm2Ag')
corrcoef(lm1',lm2')
tmmid = length(lm1);
nd = 14;
l1 = lm1(nd+1:tmmid-2*nd);
lagCC = zeros(2*nd+1,1);
for lg = -nd:nd
    l2 = lm2(nd+1+lg:tmmid-2*nd+lg);
    mid = corrcoef(l1',l2');
    lagCC(lg+nd+1) = mid(1,2);
end
txx2 = [1979.5:1./12:2021.49];
plot([-nd:nd],lagCC,'linewidth',1.5);hold on

ylabel('CC')
set(gca,'FontSize',8,'LineWidth',1,'FontName','Helvetica');
set(gca,'FontName','Helvetica','FontSize',8,'LineWidth',1.,'TickDir','in',...
    'XMinorTick','on','TickLength',[0.025 0.05],...
    'YMinorTick','on');
xlim([-nd,nd])
grid on
gtext('b Region-B SH and N^2','FontSize',8,'FontName','Helvetica','fontweight','bold')

%
axes(ha(3));
lm1Ag = mean(LTN20NCEP_AA90_94,'omitnan');
lm2Ag = mean(LH_AA90_94,'omitnan');
nd = 3;nd2 =7;
lm1mid = running_average(lm1Ag,nd,nd,nd2);
lm2mid = running_average(lm2Ag,nd,nd,nd2);
% lm1 = detrend(lm1mid(nd+1:tmYR-nd));
% lm2 = detrend(lm2mid(nd+1:tmYR-nd));
lm1 = (lm1mid(nd+1:tmYR-nd));
lm2 = (lm2mid(nd+1:tmYR-nd));
corrcoef(lm1Ag',lm2Ag')
corrcoef(lm1',lm2')
tmmid = length(lm1);
nd = 14;
l1 = lm1(nd+1:tmmid-2*nd);
lagCC = zeros(2*nd+1,1);
for lg = -nd:nd
    l2 = lm2(nd+1+lg:tmmid-2*nd+lg);
    mid = corrcoef(l1',l2');
    lagCC(lg+nd+1) = mid(1,2);
end
txx2 = [1979.5:1./12:2021.49];
plot([-nd:nd],lagCC,'linewidth',1.5);hold on

ylabel('CC')
set(gca,'FontSize',8,'LineWidth',1,'FontName','Helvetica');
set(gca,'FontName','Helvetica','FontSize',8,'LineWidth',1.,'TickDir','in',...
    'XMinorTick','on','TickLength',[0.025 0.05],...
    'YMinorTick','on');
xlim([-nd,nd])
grid on
gtext('c Region-A LH and N^2','FontSize',8,'FontName','Helvetica','fontweight','bold')

%
axes(ha(4));
lm1Ag = mean(LTN20NCEP_BB90_94,'omitnan');
lm2Ag = mean(LH_BB90_94,'omitnan');
nd = 3;nd2 =7;
lm1mid = running_average(lm1Ag,nd,nd,nd2);
lm2mid = running_average(lm2Ag,nd,nd,nd2);
% lm1 = detrend(lm1mid(nd+1:tmYR-nd));
% lm2 = detrend(lm2mid(nd+1:tmYR-nd));
lm1 = (lm1mid(nd+1:tmYR-nd));
lm2 = (lm2mid(nd+1:tmYR-nd));
corrcoef(lm1Ag',lm2Ag')
corrcoef(lm1',lm2')
tmmid = length(lm1);
nd = 14;
l1 = lm1(nd+1:tmmid-2*nd);
lagCC = zeros(2*nd+1,1);
for lg = -nd:nd
    l2 = lm2(nd+1+lg:tmmid-2*nd+lg);
    mid = corrcoef(l1',l2');
    lagCC(lg+nd+1) = mid(1,2);
end
txx2 = [1979.5:1./12:2021.49];
plot([-nd:nd],lagCC,'linewidth',1.5);hold on

ylabel('CC')
set(gca,'FontSize',8,'LineWidth',1,'FontName','Helvetica');
set(gca,'FontName','Helvetica','FontSize',8,'LineWidth',1.,'TickDir','in',...
    'XMinorTick','on','TickLength',[0.025 0.05],...
    'YMinorTick','on');
xlim([-nd,nd])
grid on
gtext('d Region-B LH and N^2','FontSize',8,'FontName','Helvetica','fontweight','bold')

axes(ha(5));
lm1Ag = mean(LTN20NCEP_AA90_94,'omitnan');
lm2Ag = mean(omega_AA90_94,'omitnan');
nd = 3;nd2 =7;
lm1mid = running_average(lm1Ag,nd,nd,nd2);
lm2mid = running_average(lm2Ag,nd,nd,nd2);
% lm1 = detrend(lm1mid(nd+1:tmYR-nd));
% lm2 = detrend(lm2mid(nd+1:tmYR-nd));
lm1 = (lm1mid(nd+1:tmYR-nd));
lm2 = (lm2mid(nd+1:tmYR-nd));
corrcoef(lm1Ag',lm2Ag')
corrcoef(lm1',lm2')
tmmid = length(lm1);
nd = 14;
l1 = lm1(nd+1:tmmid-2*nd);
lagCC = zeros(2*nd+1,1);
for lg = -nd:nd
    l2 = lm2(nd+1+lg:tmmid-2*nd+lg);
    mid = corrcoef(l1',l2');
    lagCC(lg+nd+1) = mid(1,2);
end
txx2 = [1979.5:1./12:2021.49];
plot([-nd:nd],lagCC,'linewidth',1.5);hold on
xlabel('Lag Days')
ylabel('CC')
set(gca,'FontSize',8,'LineWidth',1,'FontName','Helvetica');
set(gca,'FontName','Helvetica','FontSize',8,'LineWidth',1.,'TickDir','in',...
    'XMinorTick','on','TickLength',[0.025 0.05],...
    'YMinorTick','on');
xlim([-nd,nd])
grid on
gtext('e Region-A \omega and N^2','FontSize',8,'FontName','Helvetica','fontweight','bold')

axes(ha(6));
lm1Ag = mean(omega_BB90_94,'omitnan');
lm2Ag = mean(LTN20NCEP_BB90_94,'omitnan');
nd = 3;nd2 =7;
lm1mid = running_average(lm1Ag,nd,nd,nd2);
lm2mid = running_average(lm2Ag,nd,nd,nd2);
% lm1 = detrend(lm1mid(nd+1:tmYR-nd));
% lm2 = detrend(lm2mid(nd+1:tmYR-nd));
lm1 = (lm1mid(nd+1:tmYR-nd));
lm2 = (lm2mid(nd+1:tmYR-nd));
corrcoef(lm1Ag',lm2Ag')
corrcoef(lm1',lm2')
tmmid = length(lm1);
nd = 14;
l1 = lm1(nd+1:tmmid-2*nd);
lagCC = zeros(2*nd+1,1);
for lg = -nd:nd
    l2 = lm2(nd+1+lg:tmmid-2*nd+lg);
    mid = corrcoef(l1',l2');
    lagCC(lg+nd+1) = mid(1,2);
end
txx2 = [1979.5:1./12:2021.49];
plot([-nd:nd],lagCC,'linewidth',1.5);hold on

xlabel('Lag Days')
ylabel('CC')
set(gca,'FontSize',8,'LineWidth',1,'FontName','Helvetica');
set(gca,'FontName','Helvetica','FontSize',8,'LineWidth',1.,'TickDir','in',...
    'XMinorTick','on','TickLength',[0.025 0.05],...
    'YMinorTick','on');
xlim([-nd,nd])
grid on
gtext('f Region-B \omega and N^2','FontSize',8,'FontName','Helvetica','fontweight','bold')

sav_name=['npj/final submission/Fig7Lag3factorsRG7d71990-1994All'];
print ('-r600','-djpeg',(sav_name))
print ('-r600','-depsc2',(sav_name))
print ('-r600','-dpdf',(sav_name))
%% the table 
load('heatmonthly.mat')
NCEPDATA = load('DAY49/N2andTheta3.mat');
thetaNCEP= NCEPDATA.theta;
N20NCEP = NCEPDATA.N20;
fid = 'NCEPdata/air.mon.mean.nc';
ncdisp(fid)
airNCEP =  ncread(fid,'air');
levelNCEP = ncread(fid,'level');
[xmNCEP,ymNCEP,zmNCEP,tmNCEP] = size(airNCEP);
fid = 'NCEPdata/land.nc';
ncdisp(fid)
landNCEP = ncread(fid,'land');% 1 is land; 0 is ocean
lon0NCEP = double(ncread(fid,'lon'));
latNCEP = double(ncread(fid,'lat'));
lat = latNCEP;
lon0 = lon0NCEP;
LTN20NCEP = squeeze(mean(N20NCEP(:,:,2:3,:),3));%700-850,850-925
% === multiyear average and anomaly
[xm,ym,zm,tm] = size(omega);
mairNCEP = zeros(xmNCEP,ymNCEP,zmNCEP,12);
dairNCEP = zeros(xmNCEP,ymNCEP,zmNCEP,tmNCEP);
mLH = zeros(xm,ym,12);
dLH = zeros(xm,ym,tm);
mSH = zeros(xm,ym,12);
dSH = zeros(xm,ym,tm);
mUL = zeros(xm,ym,12);
dUL = zeros(xm,ym,tm);
mDL = zeros(xm,ym,12);
dDL = zeros(xm,ym,tm);
mUS = zeros(xm,ym,12);
dUS = zeros(xm,ym,tm);
mDS = zeros(xm,ym,12);
dDS = zeros(xm,ym,tm);
momega = zeros(xm,ym,zm,12);
domega = zeros(xm,ym,zm,tm);
diffT = squeeze(thetaNCEP(:,:,4,:) -thetaNCEP(:,:,2,:)); %850 -925 theta 
mdiffT= zeros(xm,ym,12);
ddiffT = zeros(xm,ym,tm);
mLTN20NCEP = zeros(xmNCEP,ymNCEP,12);
dLTN20NCEP = zeros(xmNCEP,ymNCEP,tmNCEP);
for jt =1:12
    mairNCEP(:,:,:,jt) = mean(airNCEP(:,:,:,jt:12:360),4);
    dairNCEP(:,:,:,jt:12:end) = airNCEP(:,:,:,jt:12:end)-mairNCEP(:,:,:,jt);
    mLH(:,:,jt) = mean(LHlrNCEP(:,:,jt:12:360),3);
    dLH(:,:,jt:12:end) = LHlrNCEP(:,:,jt:12:end)-mLH(:,:,jt);
    mSH(:,:,jt) = mean(SHlrNCEP(:,:,jt:12:360),3);
    dSH(:,:,jt:12:end) = SHlrNCEP(:,:,jt:12:end)-mSH(:,:,jt);
    mUL(:,:,jt) = mean(ULlrNCEP(:,:,jt:12:360),3);
    dUL(:,:,jt:12:end) = ULlrNCEP(:,:,jt:12:end)-mUL(:,:,jt);
    mDL(:,:,jt) = mean(DLlrNCEP(:,:,jt:12:360),3);
    dDL(:,:,jt:12:end) = DLlrNCEP(:,:,jt:12:end)-mDL(:,:,jt);
    mUS(:,:,jt) = mean(USlrNCEP(:,:,jt:12:360),3);
    dUS(:,:,jt:12:end) = USlrNCEP(:,:,jt:12:end)-mUS(:,:,jt);
    mDS(:,:,jt) = mean(DSlrNCEP(:,:,jt:12:360),3);
    dDS(:,:,jt:12:end) = DSlrNCEP(:,:,jt:12:end)-mDS(:,:,jt);
    momega(:,:,:,jt) = mean(omega(:,:,:,jt:12:360),4);
    domega(:,:,:,jt:12:end) = omega(:,:,:,jt:12:end)-momega(:,:,:,jt);
    mdiffT(:,:,jt) = mean(diffT(:,:,jt:12:360),3);
    ddiffT(:,:,jt:12:end) = diffT(:,:,jt:12:end)-mdiffT(:,:,jt);
    mLTN20NCEP(:,:,jt) = mean(LTN20NCEP(:,:,jt:12:360),3);
    dLTN20NCEP(:,:,jt:12:end) = LTN20NCEP(:,:,jt:12:end)-mLTN20NCEP(:,:,jt);
end
%% pick up the points and regional average
[xmNCEP,ymNCEP,zmNCEP,tmNCEP] = size(omega);
dLTN20NCEP_AA = ones(xmNCEP*7,tmERA5).*nan;
dLTN20NCEP_BB = ones(xmNCEP*7,tmERA5).*nan;
dair925_AA = ones(xmNCEP*7, tmERA5).*nan;
dair925_BB = ones(xmNCEP*7, tmERA5).*nan;
dair850_AA = ones(xmNCEP*7, tmERA5).*nan;
dair850_BB = ones(xmNCEP*7, tmERA5).*nan;
domega_AA = ones(xmNCEP*7,tmERA5).*nan;
domega_BB = ones(xmNCEP*7,tmERA5).*nan;
dSH_AA = ones(xmNCEP*7, tmERA5).*nan;
dSH_BB = ones(xmNCEP*7, tmERA5).*nan;
dLH_AA = ones(xmNCEP*7, tmERA5).*nan;
dLH_BB = ones(xmNCEP*7, tmERA5).*nan;
dUL_AA = ones(xmNCEP*7, tmERA5).*nan;
dDL_AA = ones(xmNCEP*7, tmERA5).*nan;
dUS_AA = ones(xmNCEP*7, tmERA5).*nan;
dDS_AA = ones(xmNCEP*7, tmERA5).*nan;
dUL_BB = ones(xmNCEP*7, tmERA5).*nan;
dDL_BB = ones(xmNCEP*7, tmERA5).*nan;
dUS_BB = ones(xmNCEP*7, tmERA5).*nan;
dDS_BB = ones(xmNCEP*7, tmERA5).*nan;
ddiffT_AA = ones(xmNCEP*7,tmERA5).*nan;
ddiffT_BB = ones(xmNCEP*7,tmERA5).*nan;
kk1 = 1;
kk2 = 1;
for i = 1:xmNCEP
    for j = 1:25
        % === pick up the points with the trend
        ln1mNCEP = (squeeze(dLTN20NCEP(i,j,1:tmERA5)));
        ln3mNCEP = (squeeze(dairNCEP(i,j,2,1:tmERA5)));
        ln4mNCEP = (squeeze(dairNCEP(i,j,3,1:tmERA5)));
        ln5mNCEP = (squeeze(domega(i,j,6,1:tmERA5)));%500hPa
        ln6mNCEP = (squeeze(dLH(i,j,1:tmERA5)));
        ln7mNCEP = (squeeze(dSH(i,j,1:tmERA5)));
        ln8mNCEP = (squeeze(dUL(i,j,1:tmERA5)));
        ln9mNCEP = (squeeze(dDL(i,j,1:tmERA5)));
        ln10mNCEP = (squeeze(dUS(i,j,1:tmERA5)));
        ln11mNCEP = (squeeze(dDS(i,j,1:tmERA5)));
        ln12mNCEP = (squeeze(ddiffT(i,j,1:tmERA5)));
        if (maskNor(i,j)==3)%norwegian
            dLTN20NCEP_AA(kk1,:) = ln1mNCEP;
            dair925_AA(kk1,:)  = ln3mNCEP;
            dair850_AA(kk1,:)  = ln4mNCEP;
            domega_AA(kk1,:)  = ln5mNCEP;
            dLH_AA(kk1,:)  = ln6mNCEP;
            dSH_AA(kk1,:)  = ln7mNCEP;
            dUL_AA(kk1,:) = ln8mNCEP;
            dDL_AA(kk1,:) = ln9mNCEP;
            dUS_AA(kk1,:) = ln10mNCEP;
            dDS_AA(kk1,:) = ln11mNCEP;
            ddiffT_AA(kk1,:)  = ln12mNCEP;
            kk1 = kk1+1;
        elseif (maskNor(i,j)==1)%iceland
            dLTN20NCEP_BB(kk2,:) = ln1mNCEP;
            dair925_BB(kk2,:) = ln3mNCEP;
            dair850_BB(kk2,:) = ln4mNCEP;
            domega_BB(kk2,:) = ln5mNCEP;
            dLH_BB(kk2,:) = ln6mNCEP;
            dSH_BB(kk2,:) = ln7mNCEP;
            dUL_BB(kk2,:) = ln8mNCEP;
            dDL_BB(kk2,:) = ln9mNCEP;
            dUS_BB(kk2,:) = ln10mNCEP;
            dDS_BB(kk2,:) = ln11mNCEP;
            ddiffT_BB(kk2,:)  = ln12mNCEP;
            kk2 = kk2+1;
        end
    end
end


%% ==== SH AND LH AND N2 REGION-A
nd = 6;nd2 =13;
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
lm1Ag = detrend(mean(dSH_AA,'omitnan'));
lm2Ag = detrend(mean(dLTN20NCEP_AA,'omitnan'));
corrcoef(lm1Ag',lm2Ag')
%% ====
nd = 6;nd2 =13;
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
lm1Ag = detrend(mean(dLH_AA,'omitnan'));
lm2Ag = detrend(mean(dLTN20NCEP_AA,'omitnan'));
corrcoef(lm1Ag',lm2Ag')

%%
nd = 6; nd2 = 13;
lm1Ag = detrend(mean(dLTN20NCEP_AA,'omitnan'));
lm2Ag = detrend(mean(domega_AA,'omitnan'));
lm1mid = running_average(lm1Ag,nd,nd,nd2);
lm2mid = running_average(lm2Ag,nd,nd,nd2);
lm1 = (lm1mid(nd+1:tmERA5-nd));
lm2 = (lm2mid(nd+1:tmERA5-nd));
corrcoef(lm1',lm2')

%% ==== SH AND LH AND N2 REGION-B
nd = 6;nd2 =13;
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
lm1Ag = detrend(mean(dSH_BB,'omitnan'));
lm2Ag = detrend(mean(dLTN20NCEP_BB,'omitnan'));
corrcoef(lm1Ag',lm2Ag')

%% ====
nd = 6;nd2 =13;
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
lm1Ag = detrend(mean(dLH_BB,'omitnan'));
lm2Ag = detrend(mean(dLTN20NCEP_BB,'omitnan'));
corrcoef(lm1Ag',lm2Ag')

%%
nd = 6; nd2 = 13;
lm1Ag = detrend(mean(dLTN20NCEP_BB,'omitnan'));
lm2Ag = detrend(mean(domega_BB,'omitnan'));
lm1mid = running_average(lm1Ag,nd,nd,nd2);
lm2mid = running_average(lm2Ag,nd,nd,nd2);
lm1 = (lm1mid(nd+1:tmERA5-nd));
lm2 = (lm2mid(nd+1:tmERA5-nd));
corrcoef(lm1',lm2')


