% Run SSW_bootstrap.m before this!
Exps={'WA1Jan35_EQ1M40_','WA1Jan35_EQ1M40lon60180_','WA1Jan35_EQ1M40lon180300_','WA1Jan35_EQ1M40lon30060_'};
lonranges=[0,360;0,360;0,360;0,360];
Amps=[0,1,2,3,4,5,6,7,10];
% Exps={'WA1Jan35_EQ1M40lon60180_','WA1Jan35_EQ1M40lon180300_','WA1Jan35_EQ1M40lon30060_'};
% Amps=[0,1,2,3,4,5,6,7,10,30];

cA={'0','1','2','3','4','5','6','7','10'};
%cA={'0','1','2','3','4','5','6','7','10','30'};
NExps=length(Exps);
colors={'k','r','b','g'};
%colors={'r','b','g'};
N=length(Amps);
nphase=40;
IDEAL=1;
boxlat0=50;
boxlat1=90;
boxlev0=600;
boxlev1=100;

% transmission ~ convergence
%  boxlat0=55;
%  boxlat1=80;
%  boxlev0=500;
%  boxlev1=80;

MJOlat=30;
calcorr=true;


%% initialization
addpath('/Users/wanyingkang/Desktop/research/matlab/code/basic/')
nSSW=zeros(NExps,N,1);
uni_power1=zeros(NExps,N,1);
notunif1=zeros(NExps,N,1);
Eady1=zeros(NExps,N,1);
U30=zeros(NExps,1);
bsyn=zeros(NExps,1);
r_bsyn=zeros(NExps,1);
U300=zeros(NExps,1);
b=zeros(NExps,1);
r_b=zeros(NExps,1);
T_phase1=zeros(NExps,N,nphase);
T_phase2=zeros(NExps,N,nphase);
U_phase1=zeros(NExps,N,nphase);
T_mean1=zeros(NExps,N,1);
U_mean1=zeros(NExps,N,1);
U_mean2=zeros(NExps,N,1);
EP=zeros(NExps,N,1);
EPstrat=zeros(NExps,N,1);
EPsta=zeros(NExps,N,1);
EPsta_ux=zeros(NExps,N,1);
EPstaint_ux=zeros(NExps,N,1);
EPstak2int_ux=zeros(NExps,N,1);
EPstak1int_ux=zeros(NExps,N,1);
EP40=zeros(NExps,N,1);
EPMJO=zeros(NExps,N,1);
EPw=zeros(NExps,N,1);
EPlk2=zeros(NExps,N,1);
EPMJOint=zeros(NExps,N,1);
EPstaint=zeros(NExps,N,1);
EPstak2=zeros(NExps,N,1);
EPstak2int=zeros(NExps,N,1);
EPstak1int=zeros(NExps,N,1);
EPstratint=zeros(NExps,N,1);
EPstratk1int=zeros(NExps,N,1);
EPstratk2int=zeros(NExps,N,1);
EPstratk2=zeros(NExps,N,1);
EPstratstak1int=zeros(NExps,N,1);
EPstratstak2int=zeros(NExps,N,1);
EPstratstak2=zeros(NExps,N,1);
EPstratstaint=zeros(NExps,N,1);
EPstratsta=zeros(NExps,N,1);
EPwint=zeros(NExps,N,1);
EP40int=zeros(NExps,N,1);
EP40int_T=zeros(NExps,N,1);
EP40int_B=zeros(NExps,N,1);
EP40int_S=zeros(NExps,N,1);
EP40int_N=zeros(NExps,N,1);
EP40stratint=zeros(NExps,N,1);
EP40strat=zeros(NExps,N,1);
EP40construint=zeros(NExps,N,1);
EP40construstratint=zeros(NExps,N,1);
EP40k2int=zeros(NExps,N,1);
EP40k2stratint=zeros(NExps,N,1);
EP40k2strat=zeros(NExps,N,1);
EP40k1int=zeros(NExps,N,1);
EP40k1stratint=zeros(NExps,N,1);
EP40k1strat=zeros(NExps,N,1);
pene=zeros(NExps,N,1);
peneint=zeros(NExps,N,1);
EPlk2=zeros(NExps,N,1);
EPlk1=zeros(NExps,N,1);
EPl=zeros(NExps,N,1);
EPint=zeros(NExps,N,1);
EPint_T=zeros(NExps,N,1);
EPint_B=zeros(NExps,N,1);
EPint_S=zeros(NExps,N,1);
EPint_N=zeros(NExps,N,1);
EPint_stratpene=zeros(NExps,N,1);
EPlint_stratpene=zeros(NExps,N,1);
EPlk2int=zeros(NExps,N,1);
EPlk1int=zeros(NExps,N,1);
EPlint=zeros(NExps,N,1);
Utropo=zeros(NExps,N,1);
Utropoint=zeros(NExps,N,1);
U_Amp_strat=zeros(NExps,N,1);
U_Amp_mid1=zeros(NExps,N,1);
U_Amp_mid2=zeros(NExps,N,1);
U_Amp_subtrop=zeros(NExps,N,1);
Cov_exit=zeros(NExps,N,1);
Uk1=zeros(NExps,N,1);
Utropolon1=zeros(NExps,N,144);
U_xavg=zeros(NExps,N,96);
jetrange=zeros(NExps,N,2);
Ks=zeros(NExps,N,1);
imaxes=zeros(NExps,N);
alphal_syn=zeros(NExps,1);
intcpt=zeros(NExps,1);
r_alphal=zeros(NExps,1);
Cu=zeros(NExps,1);
U20=zeros(NExps,1);
r_Cu=zeros(NExps,1);


%zonal asymmetry
lat_conv=40.;
lev_conv=250.;

%Ks
lev_Ks=30.;
lat_Ks=75.;

%notunif
lev1=20.;
%notunif at troposphere
lev3=50.;

%proj
lev_proj_strat=10.;
lat_proj_strat=80.;
lev_proj_mid1=50.;
lat_proj_mid1=63.;

lev_proj_subtrop=300.;
lat_proj_subtrop=30.;
lev_proj_mid2=80.;
lat_proj_mid2=50.;

%EP
levEP=boxlev1;
latEP=55.;
latEPstrat=75.;
levEPstrat=10.;

%levEP40int=lev_proj_mid2;
levEP40int=boxlev1;
latEP40int0=boxlat0;
latEP40int1=boxlat1;

levEPsta=boxlev1;
latEPsta=55.;
latEPstaint0=boxlat0;
latEPstaint1=boxlat1;
levEPsta_ux=400;
latEPstaint0_ux=40;
latEPstaint1_ux=55;
latEPsta_ux=45;
latEPstak2=45.;
levEPstak2=300.;
levEPint_stratpene=60.;

levEPstratint=levEPstrat;
latEPstratint0=50.;
latEPstratint1=90.;

levEPpene=10.;

levEPMJO=300.;
latEPMJO=MJOlat;
levEPMJOint0=boxlev0;
levEPMJOint1=boxlev1;

%Eady
lat0=50;
lev0=500;
delta_lat=6;
delta_lev=100;
lat_Ueady=75.;
lev_Ueady=300.;
lat_Ueadyint0=boxlat0;
lat_Ueadyint1=boxlat1;
lev_Ueadyint0=boxlev0;
lev_Ueadyint1=boxlev1;

a=6371000.;
deltaytropo=(latEPstaint1-latEPstaint0)/180*pi*6371000;
deltaysta=(latEPstaint1-latEPstaint0)/180*pi*6371000;
deltaystrat=(latEPstratint1-latEPstratint0)./180.*pi*6371000.;
deltap=100.*(levEPMJOint0-levEPMJOint1);

%name
cE={'A','B','C','D'};
c=strcat(repmat(cE',[1,length(cA)]),repmat(cA,[length(cE),1]));


%% SSW copy from SSW_bootstrap.m
nSSW=squeeze(mean(SSW_cis1(:,:,1:NExps),2))';

%% signal transmitted to the Arctic stratopshere
% figure4's Amp
for iExp=1:NExps
for icase=1:N
file1=char(strcat('./T_40day_',Exps(iExp),sprintf('%d',Amps(icase)),'.nc'));
U__=double(squeeze(mean(ncread(file1,'T_Amp'))))';
try
    lev=double(ncread(file1,'lev'));
end
try
    lev=double(ncread(file1,'p'));
end
lat=double(ncread(file1,'lat'));
deltay1=(lat(2)-lat(1))/180*pi*a;
[~,ilev_strat]=min(abs(lev-lev_proj_strat));
[~,ilat_strat]=min(abs(lat-lat_proj_strat));
U_Amp_strat(iExp,icase)=max(max(U__(ilev_strat-5:ilev_strat+5,ilat_strat-5:ilat_strat+5)));
[~,ilev_mid]=min(abs(lev-lev_proj_mid1));
[~,ilat_mid]=min(abs(lat-lat_proj_mid1));
U_Amp_mid1(iExp,icase)=mean(mean(U__(ilev_mid-1:ilev_mid+1,ilat_mid-2:ilat_mid+2)));
[~,ilev_mid]=min(abs(lev-lev_proj_mid2));
[~,ilat_mid]=min(abs(lat-lat_proj_mid2));
U_Amp_mid2(iExp,icase)=mean(mean(U__(ilev_mid-2:ilev_mid+2,ilat_mid-4:ilat_mid+4)));
[~,ilev_subtrop]=min(abs(lev-lev_proj_subtrop));
[~,ilat_subtrop]=min(abs(lat-lat_proj_subtrop));
U_Amp_subtrop(iExp,icase)=mean(mean(U__(ilev_subtrop-2:ilev_subtrop+2,ilat_subtrop-4:ilat_subtrop+4)));
end

%U_Amp_subtrop(iExp,:)=U_Amp_subtrop(iExp,:)-U_Amp_subtrop(iExp,1);
%U_Amp_mid1(iExp,:)=U_Amp_mid1(iExp,:)-U_Amp_mid1(iExp,1);
%U_Amp_mid2(iExp,:)=U_Amp_mid2(iExp,:)-U_Amp_mid2(iExp,1);
%U_Amp_strat(iExp,:)=U_Amp_strat(iExp,:)-U_Amp_strat(iExp,1);
end
% figure(2)
% clf
% for iExp=1:NExps
% scatter(U_Amp_strat(iExp,:),nSSW(iExp,:),char(strcat(colors(iExp),'*')));
% hold on
% text(U_Amp_strat(iExp,:)+0.001,(nSSW(iExp,:))+0.001,c(iExp,:),'FontSize',12,'Color',char(colors(iExp)))
% end
% xlabel('T Amp @ Arctic Strato.')
% ylabel('#SSW')
% set(gca,'FontSize',18)

%% waves generated in mid-lat
for iExp=1:NExps
for icase=1:N  
    file1=char(strcat('./All_',Exps(iExp),sprintf('%d',Amps(icase)),'.nc'));
    Utropo__=double(ncread(file1,'U'));
    Utropo_=squeeze(mean(Utropo__))';
    lat=double(squeeze(ncread(file1,'lat')));
    lev=double(squeeze(ncread(file1,'lev')));
    lon=double(squeeze(ncread(file1,'lon')));
    dp=double(squeeze(ncread(file1,'ilev')));
    dp=(dp(2:end)-dp(1:end-1)).*100.;
    [~,ilat]=min(abs(lat-lat_Ueady));
    [~,ilev]=min(abs(lev-lev_Ueady));
    Utropo(iExp,icase)=mean(mean(Utropo_(ilev-1:ilev+1,ilat-7:ilat+7)));
    ilat=find(lat>lat_Ueadyint0&lat<lat_Ueadyint1);
    ilev=find(lev<lev_Ueadyint0&lev>lev_Ueadyint1);
    Utropoint(iExp,icase)=sum(sum(Utropo_(ilev,ilat).*repmat(dp(ilev),[1,length(ilat)]).*repmat(cos(lat(ilat)/180*pi)',[length(ilev),1])))/sum(dp(ilev))/sum(cos(lat(ilat)/180*pi)); 
    %Utropoint(iExp,icase)=max(max(Utropo_(ilev,ilat))); 
end
end
% figure(3)
% clf
% for iExp=1:NExps
% scatter(Utropoint(iExp,:),nSSW(iExp,:),char(strcat(colors(iExp),'*')));
% hold on
% text(Utropoint(iExp,:)+0.001,(nSSW(iExp,:))+0.001,c(iExp,:),'FontSize',12,'Color',char(colors(iExp)))
% end
% xlabel('mid-lat jet strength')
% ylabel('#SSW')
% set(gca,'FontSize',18)

%% zonal asymmetry / jet exit convergence
for iExp=1:NExps
for icase=1:N  
    file1=char(strcat('./All_',Exps(iExp),sprintf('%d',Amps(icase)),'.nc'));
    Utropo__=double(ncread(file1,'U'));
    [~,ilat]=min(abs(lat-lat_conv));
    [~,ilev]=min(abs(lev-lev_conv));
    r1=-11;
    r2=11;
    %Utropolon=smoothn(squeeze(mean(Utropo__(:,ilat+r1:ilat+r2,ilev-1:ilev+1),3)),'robust');
    Utropolon=squeeze(mean(Utropo__(:,ilat+r1:ilat+r2,ilev-1:ilev+1),3));
    %Cov_exit(icase)=max(max(Utropolon(1:end-1,:)-Utropolon(2:end,:))./cos(lat(ilat+r1:ilat+r2)/180*pi)');
    %Uk1(icase)=max(abs(sum(Utropolon.*repmat(exp(i*lon/180*pi),[1,r2-r1+1]),1))/length(lon)./cos(lat(ilat+r1:ilat+r2)/180*pi)');
    if lonranges(iExp,1)<lonranges(iExp,2)
    ilon=find(lon>=lonranges(iExp,1) & lon<=lonranges(iExp,2));
    else
    ilon=[find(lon<=lonranges(iExp,1)); find(lon>=lonranges(iExp,2))];
    end
    %Utropolon_lonrange=Utropolon(ilon,:);
    [~,imax_U]=max(max(Utropolon,[],1));
    %[Cov_exit(iExp,icase),imax_conv]=max(max(Utropolon_lonrange(1:end-1,:)-Utropolon_lonrange(2:end,:)));
    [Cov_exit(iExp,icase),imax_conv]=max(max(Utropolon(1:end-1,:)-Utropolon(2:end,:)));
    %Utropolon1(iExp,icase,:)=Utropolon(:,imax_conv);
    %Utropolon1(iExp,icase,:)=Utropolon(:,imax_U);
    [Utropolon1(iExp,icase,:),ilat_maxU]=max(Utropolon,[],2);
    lat1=lat(ilat+r1:ilat+r2);
    lonband1=200;
    lonband0=60;
    ilonband=(lon>lonband0&lon<lonband1);
    jetrange(iExp,icase,:)=[lat1(min(ilat_maxU(ilonband))),lat1(max(ilat_maxU(ilonband)))];
    %Utropolon1(iExp,icase,:)=Utropolon(:,1);
    imaxes(iExp,icase)=imax_conv;
    
    %Utropolon1(iExp,icase,:)=squeeze(mean(Utropolon,2));
    Uk1(iExp,icase)=mean(abs(sum(Utropolon.*repmat(exp(1i*lon/180*pi),[1,r2-r1+1]),1))/length(lon));
    
    %U_xavg
    [~,ilev]=min(abs(lev-lev_conv));
    U_xavg(iExp,icase,:)=squeeze(mean(mean(Utropo__(:,:,ilev-1:ilev+1)),3));
end
end

% figure(16)
% clf
% for iExp=1:NExps
% scatter(Uk1(iExp,:),EPstaint_ux(iExp,:),char(strcat(colors(iExp),'*')))
% hold on
% end
% xlabel(strcat('U k=1 zonal asymmetry @',sprintf('%d',lat_conv),'N'))
% ylabel(strcat('EP stationary ',sprintf('%d',latEPstaint0_ux),'-',sprintf('%d',latEPstaint1_ux),'N'))
% set(gca,'FontSize',18)

disp('k=1 U ampitude seems to correlate with the upward stationary EP in the same region  (EPstaint_ux)')
disp('Uk1 ~ EPstaint also seems acceptable')


% figure(4)
% clf
% for iExp=1:NExps
% scatter(Uk1(iExp,:),nSSW(iExp,:),char(strcat(colors(iExp),'*')));
% hold on
% text(Uk1(iExp,:)+0.001,(nSSW(iExp,:))+0.001,c(iExp,:),'FontSize',12,'Color',char(colors(iExp)))
% end
% xlabel('U wind convergence')
% ylabel('#SSW')
% set(gca,'FontSize',18)



if calcorr
%% calculate EPs
for iExp=1:NExps
for icase=1:N  
    file1=char(strcat('./yz_EP_filter_010000_DJF_',Exps(iExp),sprintf('%d',Amps(icase)),'k2.nc'));
    EP_=-double(squeeze(ncread(file1,'F_k')))';
    lat=double(squeeze(ncread(file1,'lat')));
    lev=double(squeeze(ncread(file1,'lev')));
    [~,ilat]=min(abs(lat-latEPsta));
    [~,ilev]=min(abs(lev-levEPsta));
    EPlk2(iExp,icase)=mean(mean(EP_(ilev-1:ilev+1,ilat-1:ilat+1))).*cos(lat(ilat)*pi/180.); 
    ilat=find((lat>latEPstaint0& lat<latEPstaint1));
    EPlk2int(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    EPlk2int(iExp,icase)=EPlk2int(iExp,icase); %stationary wave included
    [~,ilev]=min(abs(lev-levEPstrat));
    ilat=find((lat>latEPstratint0& lat<latEPstratint1));
    EPstratk2int(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    EPstratk2int(iExp,icase)=EPstratk2int(iExp,icase);
    [~,ilat]=min(abs(lat-latEPstrat));
    EPstratk2(iExp,icase)=mean(mean(EP_(ilev-1:ilev+1,ilat-1:ilat+1))).*cos(lat(ilat)*pi/180.); 
    
    file1=char(strcat('./yz_EP_filter_010000_DJF_',Exps(iExp),sprintf('%d',Amps(icase)),'k1.nc'));
    EP_=-double(squeeze(ncread(file1,'F_k')))';
    lat=double(squeeze(ncread(file1,'lat')));
    lev=double(squeeze(ncread(file1,'lev')));
    [~,ilat]=min(abs(lat-latEPsta));
    [~,ilev]=min(abs(lev-levEPsta));
    EPlk1(iExp,icase)=mean(mean(EP_(ilev-1:ilev+1,ilat-1:ilat+1))); 
    ilat=find((lat>latEPstaint0& lat<latEPstaint1));
    EPlk1int(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    EPlk1int(iExp,icase)=EPlk1int(iExp,icase);
    [~,ilev]=min(abs(lev-levEPstrat));
    ilat=find((lat>latEPstratint0& lat<latEPstratint1));
    EPstratk1int(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    EPstratk1int(iExp,icase)=EPstratk1int(iExp,icase);

    file1=char(strcat('./yz_EP___DJF_',Exps(iExp),sprintf('%d',Amps(icase)),'.nc'));
    EP_=-double(squeeze(ncread(file1,'F_k')))';
    EP_j=double(squeeze(ncread(file1,'F_j')))';
    lat=double(squeeze(ncread(file1,'lat')));
    lev=double(squeeze(ncread(file1,'lev')));
    ilat=find((lat>latEPstaint0& lat<latEPstaint1));
    [~,ilev]=min(abs(lev-levEPsta));
    EPint(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    
    ilat=find((lat>boxlat0& lat<boxlat1));
    [~,ilev]=min(abs(lev-boxlev1));
    EPint_T(iExp,icase)=sum(squeeze(EP_(ilev,ilat)).*deltay1);
    [~,ilev]=min(abs(lev-boxlev0));
    EPint_B(iExp,icase)=sum(squeeze(EP_(ilev,ilat)).*deltay1);
    ilev=find(lev>boxlev1&lev<boxlev0);
    [~,ilat]=min(abs(lat-boxlat0));
    EPint_S(iExp,icase)=sum(squeeze(EP_j(ilev,ilat)).*dp(ilev));
    [~,ilat]=min(abs(lat-boxlat1));
    EPint_N(iExp,icase)=sum(squeeze(EP_j(ilev,ilat)).*dp(ilev));
    
    [~,ilev]=min(abs(lev-levEPstrat));
    [~,ilat]=min(abs(lat-latEPstrat));
    EPstrat(iExp,icase)=mean(mean(EP_(ilev-1:ilev+1,ilat-1:ilat+1))).*cos(lat(ilat)*pi/180.); 

end
end
for iExp=1:NExps
for icase=1:N  
    %file1=char(strcat('./yz_EP_filter_3545_DJF_',Exps(iExp),sprintf('%d',Amps(icase)),'.nc'));
    file1=char(strcat('./yz_EP_filter_4040_DJF_',Exps(iExp),sprintf('%d',Amps(icase)),'.nc'));
    EP_j=double(squeeze(ncread(file1,'F_j')))';
    lat=double(squeeze(ncread(file1,'lat')));
    lev=double(squeeze(ncread(file1,'lev')));
    [~,ilat]=min(abs(lat-latEPMJO));
    [~,ilev]=min(abs(lev-levEPMJO));
    EPMJO(iExp,icase)=mean(mean(EP_j(ilev-1:ilev+1,ilat-2:ilat+2)));
    ilev=find((lev>levEPMJOint1& lev<levEPMJOint0));
    EPMJOint(iExp,icase)=sum(squeeze(mean(EP_j(ilev,ilat-2:ilat+2),2)).*dp(ilev)).*cos(lat(ilat)*pi/180.);
    EP_=-double(squeeze(ncread(file1,'F_k')))';
    [~,ilat]=min(abs(lat-latEPsta));
    [~,ilev]=min(abs(lev-levEPpene));
    EP40(iExp,icase)=mean(mean(EP_(ilev-2:ilev+2,ilat-2:ilat+2))).*cos(lat(ilat)*pi/180.);
    ilat=find(lat>latEP40int0&lat<latEP40int1);
    [~,ilev]=min(abs(lev-levEP40int));
    EP40int(iExp,icase)=sum(mean(EP_(ilev-2:ilev+2,ilat)).*deltay1.*cos(lat(ilat)*pi/180.)');
    
    EP40int_T(iExp,icase)=sum(squeeze(EP_(ilev,ilat)).*deltay1);
    [~,ilev]=min(abs(lev-boxlev0));
    EP40int_B(iExp,icase)=sum(squeeze(EP_(ilev,ilat)).*deltay1);
    ilev=find(lev>boxlev1&lev<boxlev0);
    [~,ilat]=min(abs(lat-boxlat0));
    EP40int_S(iExp,icase)=sum(squeeze(EP_j(ilev,ilat)).*dp(ilev));
    [~,ilat]=min(abs(lat-boxlat1));
    EP40int_N(iExp,icase)=sum(squeeze(EP_j(ilev,ilat)).*dp(ilev));

    [~,ilev]=min(abs(lev-levEPstrat));
    ilat=find((lat>latEPstratint0& lat<latEPstratint1));
    EP40stratint(iExp,icase)=sum(squeeze(mean(EP_(ilev-1:ilev+1,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    [~,ilat]=min(abs(lat-latEPstrat));
    EP40strat(iExp,icase)=mean(mean(EP_(ilev-1:ilev+1,ilat-1:ilat+1))).*cos(lat(ilat)*pi/180.); 

    
    file1=char(strcat('./yz_EP_stationary__DJF_',Exps(iExp),sprintf('%d',Amps(icase)),'.nc'));
    EP_=-double(squeeze(ncread(file1,'F_k')))';
    EP_j=double(squeeze(ncread(file1,'F_j')))';
    lat=double(squeeze(ncread(file1,'lat')));
    lev=double(squeeze(ncread(file1,'lev')));
    [~,ilat]=min(abs(lat-latEPsta));
    [~,ilev]=min(abs(lev-levEPsta));
    EPsta(iExp,icase)=mean(mean(EP_(ilev-2:ilev-+2,ilat-2:ilat+2)));
    [~,ilat]=min(abs(lat-latEPsta_ux));
    [~,ilev]=min(abs(lev-levEPsta_ux));
    r1=-5;
    r2=5;
    EPsta_ux(iExp,icase)=mean(mean(EP_(ilev-2:ilev+2,ilat+r1:ilat+r2)).*cos(lat(ilat+r1:ilat+r2).*pi./180.)');
    ilat=find((lat>latEPstaint0_ux& lat<latEPstaint1_ux));
    EPstaint_ux(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat).*pi./180.)');
    [~,ilev]=min(abs(lev-levEPsta));
    ilat=find((lat>latEPstaint0& lat<latEPstaint1));
    EPstaint(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat).*pi./180.)');
    [~,ilev]=min(abs(lev-levEPstrat));
    ilat=find((lat>latEPstratint0& lat<latEPstratint1));
    EPstratstaint(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat).*pi./180.)');
    [~,ilat]=min(abs(lat-latEPstrat));
    EPstratsta(iExp,icase)=mean(mean(EP_(ilev-1:ilev+1,ilat-1:ilat+1)).*cos(lat(ilat-1:ilat+1).*pi./180.)');
end
end

for iExp=1:NExps
for icase=1:N  
    file1=char(strcat('./yz_EP_filter_17_DJF_',Exps(iExp),sprintf('%d',Amps(icase)),'.nc'));
    EP_=-double(squeeze(ncread(file1,'F_k')))';
    lat=double(squeeze(ncread(file1,'lat')));
    lev=double(squeeze(ncread(file1,'lev')));
    [~,ilat]=min(abs(lat-latEPsta));
    [~,ilev]=min(abs(lev-levEPsta));
    EPw(iExp,icase)=mean(mean(EP_(ilev-1:ilev+1,ilat-1:ilat+1))).*cos(lat(ilat)*pi/180.)';
    ilat=find((lat>latEPstaint0& lat<latEPstaint1));
    [~,ilev]=min(abs(lev-levEPsta));
    EPwint(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
end
end
for iExp=1:NExps
for icase=1:N 
    file1=char(strcat('./yz_EP_filter_4040_DJF_',Exps(iExp),sprintf('%d',Amps(icase)),'stationary_k1.nc'));
    EP_=-double(squeeze(ncread(file1,'F_k')))';
    [~,ilev]=min(abs(lev-levEPsta));
    ilat=find((lat>latEPstaint0& lat<latEPstaint1));
    EPstak1int(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    [~,ilev]=min(abs(lev-levEPstrat));
    ilat=find((lat>latEPstratint0& lat<latEPstratint1));
    EPstratstak1int(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    [~,ilev]=min(abs(lev-levEPsta_ux));
    ilat=find((lat>latEPstaint0_ux& lat<latEPstaint1_ux));
    EPstak1int_ux(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    
    file1=char(strcat('./yz_EP_filter_4040_DJF_',Exps(iExp),sprintf('%d',Amps(icase)),'stationary_k2.nc'));
    EP_=-double(squeeze(ncread(file1,'F_k')))';
    [~,ilev]=min(abs(lev-levEPsta));
    ilat=find((lat>latEPstaint0& lat<latEPstaint1));
    EPstak2int(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    [~,ilev]=min(abs(lev-levEPstrat));
    ilat=find((lat>latEPstratint0& lat<latEPstratint1));
    EPstratstak2int(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    [~,ilat]=min(abs(lat-latEPstrat));
    EPstratstak2(iExp,icase)=mean(mean(EP_(ilev-1:ilev+1,ilat-1:ilat+1))).*cos(lat(ilat)*pi/180.); 
       
    file1=char(strcat('./yz_EP_filter_4040_DJF_',Exps(iExp),sprintf('%d',Amps(icase)),'k2.nc'));
    EP_=-double(squeeze(ncread(file1,'F_k')))';
    [~,ilev]=min(abs(lev-levEPsta));
    ilat=find((lat>latEPstaint0& lat<latEPstaint1));
    EP40k2int(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    [~,ilev]=min(abs(lev-levEPstrat));
    ilat=find((lat>latEPstratint0& lat<latEPstratint1));
    EP40k2stratint(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    [~,ilat]=min(abs(lat-latEPstrat));
    EP40k2strat(iExp,icase)=mean(mean(EP_(ilev-1:ilev+1,ilat-1:ilat+1))).*cos(lat(ilat)*pi/180.);
    
    file1=char(strcat('./yz_EP_filter_4040_DJF_',Exps(iExp),sprintf('%d',Amps(icase)),'k1.nc'));
    EP_=-double(squeeze(ncread(file1,'F_k')))';
    [~,ilev]=min(abs(lev-levEPsta));
    ilat=find((lat>latEPstaint0& lat<latEPstaint1));
    EP40k1int(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    [~,ilev]=min(abs(lev-levEPstrat));
    ilat=find((lat>latEPstratint0& lat<latEPstratint1));
    EP40k1stratint(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
    [~,ilat]=min(abs(lat-latEPstrat));
    EP40k1strat(iExp,icase)=mean(mean(EP_(ilev-1:ilev+1,ilat-1:ilat+1))).*cos(lat(ilat)*pi/180.);
    
    [~,ilev]=min(abs(lev-levEPsta_ux));
    ilat=find((lat>latEPstaint0_ux& lat<latEPstaint1_ux));
    EPstak2int_ux(iExp,icase)=sum(squeeze(mean(EP_(ilev-2:ilev+2,ilat),1)).*deltay1.*cos(lat(ilat)*pi/180.)');
end
end

%% SSW index pick (EP ~ #SSW)
EPstratfinal=EPstratk2int;
EPstratfinal=EPstratk1int;
%EPstratfinal=EPstratk2int;
EPtropofinal=EPlk2int;
EPtropofinal=EPlk1int;
%EPtropofinal=EPint;
%EPtropofinal=EPlk1int;

for iExp=1:NExps
for icase=1:N  
%constructive effect
    file1=char(strcat('./yz_EP_filter_4040_DJF_',Exps(iExp),sprintf('%d',Amps(icase)),'construct_k2.nc'));
    EP_=-double(squeeze(ncread(file1,'F_k')))';
    lat=double(squeeze(ncread(file1,'lat')));
    lev=double(squeeze(ncread(file1,'lev')));
    ilat=find(lat>latEP40int0&lat<latEP40int1);
    [~,ilev]=min(abs(lev-levEP40int));
    EP40construint(iExp,icase)=sum(mean(EP_(ilev-2:ilev+2,ilat)).*deltay1.*cos(lat(ilat)*pi/180.)');
    [~,ilev]=min(abs(lev-levEPstrat));
    ilat=find((lat>latEPstratint0& lat<latEPstratint1));
    EP40construstratint(iExp,icase)=sum(mean(EP_(ilev-1:ilev+1,ilat)).*deltay1.*cos(lat(ilat)*pi/180.)');
end

% stratospheric maximum upward EP represent #SSW
figure(30)
if iExp==1
  clf
end
plot((EPstratfinal(iExp,:)+EP40construstratint(iExp,:))/deltaystrat,nSSW(iExp,:),char(strcat(colors(iExp),'*')))
hold on
c1=c(iExp,:);
text(squeeze((EPstratfinal(iExp,:)+EP40construstratint(iExp,:))/deltaystrat),squeeze(nSSW(iExp,:)),c1,'FontSize',16,'Color',char(colors(iExp)))
xlabel('time-peak of stratospheric $EP_z$','Interpreter','latex')
ylabel('SSW number per decade')
set(gca,'FontSize',18)

% tropospheric maximum upward EP (large scale) represent #SSW
figure(31)
if iExp==1
  clf
end
plot((EPtropofinal(iExp,:)+EP40construint(iExp,:))/deltaytropo,nSSW(iExp,:),char(strcat(colors(iExp),'*')))
hold on
c1=c(iExp,:);
text((EPtropofinal(iExp,:)+EP40construint(iExp,:))/deltaytropo,nSSW(iExp,:),c1,'FontSize',16,'Color',char(colors(iExp)))
xlabel('time-peak of large scale tropospheric $EP_z$','Interpreter','latex')
ylabel('SSW number per decade')
set(gca,'FontSize',18)

% stratospheric mean upward EP represent #SSW
figure(32)
if iExp==1
  clf
end
plot((EPstratfinal(iExp,:))/deltaystrat,nSSW(iExp,:),char(strcat(colors(iExp),'*')))
hold on
c1=c(iExp,:);
text((EPstratfinal(iExp,:))/deltaystrat,nSSW(iExp,:),c1,'FontSize',16,'Color',char(colors(iExp)))
xlabel('time-mean of stratospheric $EP_z$','Interpreter','latex')
ylabel('SSW number per decade')
set(gca,'FontSize',18)

% tropospheric mean upward EP (large scale) represent #SSW
figure(33)
if iExp==1
  clf
end
plot((EPtropofinal(iExp,:))/deltaytropo,nSSW(iExp,:),char(strcat(colors(iExp),'*')))
hold on
c1=c(iExp,:);
text((EPtropofinal(iExp,:))/deltaytropo,nSSW(iExp,:),c1,'FontSize',16,'Color',char(colors(iExp)))
xlabel('time-mean of large scale tropospheric $EP_z$','Interpreter','latex')
ylabel('SSW number per decade')
set(gca,'FontSize',18)

% % tropospheric mean upward EP (all scale) represent #SSW
% figure(34)
% if iExp==1
%   clf
% end
% plot((EPint(iExp,:))/deltaytropo,nSSW(iExp,:),char(strcat(colors(iExp),'*')))
% hold on
% c1=c(iExp,:);
% text((EPint(iExp,:))/deltaytropo+dx,nSSW(iExp,:)+dy,c1,'FontSize',16,'Color',char(colors(iExp)))
% xlabel('time-mean of all scale tropospheric $EP_z$','Interpreter','latex')
% ylabel('SSW number per decade')
% set(gca,'FontSize',18)
end
disp('=======================WHICH EP FLUX USED TO PREDICT SSW?=======================')
disp('Seems the tropospheric EP is enough.')
disp('Better correlation when move boxlat1 northward (80N).')
disp('Using time-mean EP cause huge bias when MJO amplitude goes to 30K/day')
disp('Taking the constructive effect makes the 30K simulation merge with others')
disp('When using wind drop to identify SSWs, maximum EP is better index than time-mean EP.')
disp('When using wind reversal to identify SSWs, time-mean EP is better than maximum EP.')
disp('================================================================================')
disp(' ')
pause  

%% what does this time-mean upward EP in troposphere made of
for iExp=1:NExps
figure(40+iExp)
clf
nSSW(iExp,1)=nSSW(1,1);
plot(Amps, nSSW(iExp,:),'k','LineWidth',4);
ylim([-15,20])
xlim([0,10])
ylabel('SSW occurrence per decade')
ax1=gca;
ax1.YColor='k';
%set(gcf,'Position',[114   249   564   320]);
set(gcf,'Position',[114   249   400   520]);
ax1.Position=[0.13 0.15 0.7 0.75];
%[0.13 0.15 0.7 0.8];
hold on
set(gca,'FontSize',16)
ax2=axes('Position',ax1.Position,'YAxisLocation','right','Color','none','YColor','k');
%xlabel('Stationary wave at mid-lat')
hold on
showEPstrat=false;
if showEPstrat
% EPint2=EPstratk2int(iExp,:)./deltaystrat;
% EP40int2=EP40stratint(iExp,:)./deltaystrat;
% EPstaint2=EPstratstak2int(iExp,:)./deltaystrat;
EPint2=EPstrat(iExp,:)./deltaystrat;
EP40int2=EP40strat(iExp,:)./deltaystrat;
EPstaint2=EPstratsta(iExp,:)./deltaystrat;
EPconstru2=EP40construstratint(iExp,:)./deltaystrat;
else
EPint2=EPlk2int(iExp,:)./deltaytropo;
%EP40int2=EP40int(iExp,:)./deltaytropo;
EPstaint2=EPstak2int(iExp,:)./deltaytropo;
EPconstru2=EP40construint(iExp,:)./deltaytropo;
% EPint2=EPint(iExp,:)./deltaytropo;
 EP40int2=EP40k2int(iExp,:)./deltaytropo;
% EPstaint2=EPstaint(iExp,:)./deltaytropo;
end

% plot \delta EP as function of Amp.
plot(Amps, EP40int2-EP40int2(1),'r--*','LineWidth',2,'Parent',ax2);
plot(Amps, EPstaint2-EPstaint2(1),'g--*','LineWidth',2,'Parent',ax2);
plot(Amps, (EPint2-EP40int2-EPstaint2)-(EPint2(1)-EP40int2(1)-EPstaint2(1)),'b--*','LineWidth',2,'Parent',ax2);
plot(Amps, (EPint2)-(EPint2(1)),'k--*','LineWidth',2,'Parent',ax2);
lg=legend({'MJO-forced waves','stationary waves','other waves','together'});

% plot(Amps, EP40int2-EP40int2(1),'r--*','LineWidth',2,'Parent',ax2);
% plot(Amps, (EPint2-EP40int2)-(EPint2(1)-EP40int2(1)),'b--*','LineWidth',2,'Parent',ax2);
% plot(Amps, (EPint2)-(EPint2(1)),'k--*','LineWidth',2,'Parent',ax2);
% %plot(Amps, (EPconstru2)-(EPconstru2(1)),'k--o','LineWidth',1,'Parent',ax2);
% lg=legend({'MJO-forced waves','other waves','total'});



lg.Color=[1 1 1];
lg.Location='southwest';
xlim([0,10])
ylim([-0.008,0.008])
yticks([-0.008:0.002:0.008])
%plot(Amps, (EP40int+EP40construint*0.)./deltaytropo-(EP40int(1)+EP40construint(1)*0.)./deltaytropo,'r--*','LineWidth',2,'Parent',ax2);
%plot(Amps, (EPint-EPwint-EP40int)./deltaytropo-(EPint(1)-EPwint(1)-EP40int(1))./deltaytropo,'b--*','LineWidth',2,'Parent',ax2);
%ylim([-0.03,0.03])
plot([0,10],[0,0],'k--')
xlabel('MJO Amp (K/day)')
ylabel('$\Delta EP_z \ \ (m\cdot Pa/s^2)$','Interpreter','latex','Color','k')
title(char(strcat(Exps(iExp),'[Amp]')),'Interpreter','none')
set(gca,'FontSize',16)
end
disp('=======================CONTRIBUTION FROM MJO-FORCED WAVE + OTHER WAVES =======================')
disp('The 60-180E experiments actually take the least advantage of having the MJO-forced wave train transmitted.')
disp('Instead 60-180 10K experiment enhances the SSWs by enhancing the stationary waves')
disp('But how the stationary wave get strengthened by transient eddies from tropics is a question to ask!!!!')
disp('We would need the spectral model to explore the mechanism for this. The explanation in our first paper is not enough.')
disp('If using wind drop to identify SSWs, the EPflux here seems not explaining explaining the #SSW change.')
disp('===============================================================================================')
disp('  ')
pause
    
%% transmission coef.: EPsta~notunif2 / peneint
disp('=======================TRANSMISSION COEF ~ CONVERGENCE=====================')
for iExp=1:NExps
pene(iExp,:)=(EP40(iExp,:)-EP40(iExp,1)*0.0)./(EPMJO(iExp,:)-EPMJO(iExp,1)*1);
peneint(iExp,:)=((EP40int(iExp,:)-EP40int(iExp,1).*1.0))./(EPMJOint(iExp,:)-EPMJOint(iExp,1).*1.0);
%peneint(iExp,:)=EP40int(iExp,:)./EPMJOint(iExp,:);
end

figure(5)
clf
for iExp=1:NExps
%x=squeeze(EPstaint(iExp,2:end))./deltaytropo;
%x=squeeze(EPstaint_ux(iExp,2:end))./deltaytropo;
%x=squeeze(Uk1(iExp,2:end));
x=squeeze(Cov_exit(iExp,2:end));
%x=squeeze(Utropoint(iExp,2:end));
y=squeeze(peneint(iExp,2:end));
c1=c(iExp,2:end);
plot(x,y,char(strcat(colors(iExp),'*')))
hold on
dx=0.000000002;
dy=0.000001;
text(x+dx,y+dy,c1,'FontSize',16,'Color',char(colors(iExp)))

x0=mean(x);
y0=mean(y);
xa=x-x0;
ya=y-y0;
gamma02=sum(xa.*ya)/sum(xa.*xa);
Ast02=(-y0+x0*gamma02)./gamma02;
r_gamma02=min(min(corrcoef(x,y)));
hold on
xt=[min(x)-0.1*(max(x)-min(x)):0.01*(max(x)-min(x)):max(x)+0.1*(max(x)-min(x))];
yt=gamma02.*(xt-Ast02);
plot(xt,yt,char(strcat(colors(iExp),'--')))
disp(char(strcat(Exps(iExp),': gamma02 = ',sprintf('%4.3f',gamma02),', Ast02 = ',sprintf('%4.4f',Ast02),'     r = ',sprintf('%3.2f',r_gamma02))))
end
title('')
%xlabel(strcat('Stationary Eddy Upward EP Flux (',sprintf('%d',latEPsta_ux-10),'-',sprintf('%d',latEPsta_ux+10),'N, ',sprintf('%d',levEPsta_ux),'mb)'))
xlabel('Convergence at the jet exit region')
%xlabel('k=1 Amplitude of climatological U')
%xlabel('$stationary\ EP_{z} \hspace{0.5cm} (45-65N)$','Interpreter','latex')
%xlabel(strcat('Large scale Eddy Upward EP Flux (',sprintf('%d',latEPsta-10),'-',sprintf('%d',latEPsta+10),'N, ',sprintf('%d',levEPsta),'mb)'))
ylabel(strcat('$\left. \int_{',sprintf('%d',latEP40int0),'N}^{',sprintf('%d',latEP40int1),'N} \left.\it{EP}_{\it{MJO},z} \right|_{',sprintf('%d',levEP40int),'mb} dy \right/ \int_{',sprintf('%d',levEPMJOint1),'mb}^{',sprintf('%d',levEPMJOint0),'mb} \left.\it{EP}_{\it{MJO},y}\right|_{',sprintf('%d',latEPMJO),'N} dp$'),'Interpreter','latex')
set(gca,'FontSize',18)
disp('Correlation here can be improved if lift the bottom of the box to 500mb or 400mb')
disp(' ')
disp('The correlation between the jet exit convergence and transmission coef collapse well.')
disp('But the jet exit convergence is not correlated well with the averaged upward propagating waves (whatever box) or the k=1 U amplitude.')
disp('Thus we need to choose other factors than EPsta to describe the convergence.')

disp(' ')
disp('When the forcing is limited in some longitudinal range and is very strong, the zonal asymmetry will develop in one way or another somehow.')
disp('Thus the stationary wave (included in other waves) will finally increase instead of decrease as in EQ1M40s.\n')
disp('')
disp('However the convergence strength is always reduced when Amp gets stronger,')
disp('because the MJO-forced waves extract energy from the mean flow, especially 60-180E.')
disp('')

disp('============================================================================')
disp(' ')
figure(10)
clf
for iExp=1:NExps
scatter(Cov_exit(iExp,:),Utropoint(iExp,:),char(strcat(colors(iExp),'*')))
hold on
c1=c(iExp,1:end);
text(Cov_exit(iExp,:)+dx,Utropoint(iExp,:)+dy,c1,'FontSize',16,'Color',char(colors(iExp)))
end
title('Relate jet convergence to jet strength')
xlabel('Max convergence of jet')
ylabel('Mean jet speed')
set(gca,'FontSize',18)

pause

%% How efficiency MJO forcing triggers waves
figure(7)
clf
for iExp=1:NExps
x=Amps;
y=EPMJOint(iExp,:)./deltap;
plot(x,y,char(strcat(colors(iExp),'^')))
hold on
dx=0.00002;
dy=0.0001;
c1=c(iExp,:);
text(x+dx,y+dy,c1,'FontSize',16,'Color',char(colors(iExp)))
ylim([0,7])
ylabel(strcat('EP_y in MJO frequency @ ',sprintf('%d',latEPMJO),'N'))
xlabel('MJO forcing Amplitude')
set(gca,'FontSize',18)
end
disp('=======================MJO wave train strength ~ MJO Amplitude added=======================')
disp(strcat('The MJO propagation to ',sprintf('%d',latEPMJO),'N is almost linear.'))
disp('The three longitudinal restricted experiments sum up to be similar to circumglobal ones')
disp('Total waves triggered by MJO in different locations seem the same.')
disp('But different locations tend to direct the waves to upper/lower troposphere.')

disp('============================================================================================')
disp(' ')
pause

%% how zonal asymmetry is modified by waves from different longitudes
cs=jet(length(Amps));
pltdiff=true;
for iExp=1:NExps
    figure(10+iExp)
    clf
    hold on
    for icase=1:N
        yyaxis left
            set(gca,'YColor','k')
            plot(lon,squeeze(Utropolon1(iExp,icase,:)-Utropolon1(iExp,1,:)),'LineWidth',2,'Color',cs(icase,:),'LineStyle','-','Marker','none')
            ylim([-25,15])
            %ylabel(strcat('dU @',sprintf('%d',lat_conv),'N'))
            ylabel('\Delta max jet speed (m/s)')
        yyaxis right
            set(gca,'YColor','k')
            plot(lon,smooth(squeeze(Utropolon1(iExp,icase,:)),10,'lowess'),'LineWidth',2,'Color',cs(icase,:),'LineStyle','--','Marker','none')
            ylim([0,60])
            %ylabel(strcat('U @',sprintf('%d',lat_conv),'N'))
            ylabel('max jet speed (m/s)')
            jetrange1=squeeze(jetrange(iExp,icase,:));
            errorbar(360+icase*5,mean(jetrange1),(jetrange1(2)-jetrange(1))/2,'LineWidth',2,'Color',cs(icase,:));
    end
    xlabel('longitude (^\circE)')
    xticks([0:60:360])
    xlim([0,360+(N+1)*5]);
    title(strcat(char(Exps(iExp)),'Amps'),'Interpreter','none')
    set(gca,'FontSize',18)
    set(gcf,'Position',[684   517   601   369])
   
end

for iExp=1:NExps
    figure(15+iExp)
    clf
    hold on
    for icase=3:6
        yyaxis left
            plot(lon,squeeze(Utropolon1(iExp,icase,:)-Utropolon1(iExp,1,:)),'LineWidth',2,'Color',cs(icase,:),'LineStyle','--','Marker','none')
            %ylim([-25,15])
            %plot(lon,smooth(squeeze(Utropolon1(iExp,icase,:)-Utropolon1(iExp,1,:)).*squeeze(Utropolon1(iExp,icase,:)),10)./deltay1.*86400,'LineWidth',2,'Color',cs(icase,:),'LineStyle','-')
            %plot(lon,gradient(smooth(squeeze(Utropolon1(iExp,icase,:)-Utropolon1(iExp,1,:)).*squeeze(Utropolon1(iExp,icase,:)),10))./deltay1.*86400,'LineWidth',2,'Color',cs(icase,:),'LineStyle','-')
            plot(lon(2:end),diff(smooth(squeeze(Utropolon1(iExp,icase,:)-Utropolon1(iExp,1,:)).*squeeze(Utropolon1(iExp,icase,:))+0.5.*squeeze(Utropolon1(iExp,icase,:)-Utropolon1(iExp,1,:)).^2,20))./deltay1.*86400+1./80.*squeeze(Utropolon1(iExp,icase,2:end)-Utropolon1(iExp,1,2:end)),'LineWidth',2,'Color',cs(icase,:),'LineStyle','-','Marker','none')
            ylim([-3,3])
            ylabel(strcat('dU @',sprintf('%d',lat_conv),'N'))
        yyaxis right
            %plot(lon,smooth(squeeze(Utropolon1(iExp,icase,:)),10,'lowess'),'LineWidth',2,'Color',cs(icase,:),'LineStyle','--')
            %ylim([0,50])
            plot(lon,smooth(squeeze(Utropolon1(iExp,1,:)),5,'lowess'),'LineWidth',2,'Color',cs(9,:),'LineStyle','--')
            plot(lon,smooth(squeeze(Utropolon1(iExp,1,:)),50,'lowess'),'LineWidth',2,'Color',cs(1,:),'LineStyle','--')
            plot(lon,smooth(squeeze(Utropolon1(iExp,1,:)),50,'lowess')-smooth(squeeze(Utropolon1(iExp,1,:)),5,'lowess'),'LineWidth',2,'Color',cs(4,:),'LineStyle','--')
            %plot(lon,smooth(gradient(smooth(gradient(smooth(squeeze(Utropolon1(iExp,1,:)),10,'lowess')),10,'lowess')),7,'lowess'),'LineWidth',2,'Color',cs(1,:),'LineStyle','--')
            %plot(lon,squeeze(Utropolon1(iExp,1,:)),'LineWidth',2,'Color',cs(1,:),'LineStyle','--')
            %ylim([-0.3,0.3])
            ylabel(strcat('U @',sprintf('%d',lat_conv),'N'))
    end
    xlabel('longitude (^oE)')
    
    title(strcat(char(Exps(iExp)),'Amps'),'Interpreter','none')
    set(gca,'FontSize',18)
end


disp('=======================How Zonal convergence is modified by forcing=======================')
disp('It seems to show that the longitudinal range where mid-lat jet gets weakened')
disp('coincides with the forcing location.')
disp('============================================================================================')
disp(' ')
pause

%% wave generation rate, b: U2 ~ EPsyn
disp('=======================MIDLAT-GENERATED SYNOPTIC EDDY ~ JET SPEED=======================')

figure(6)
clf
for iExp=1:NExps
x=Utropo(iExp,:);
y=EPwint(iExp,:)./deltaysta;
plot(x,y,char(strcat(colors(iExp),'^')))
hold on
dx=0.00002;
dy=0.0001;
c1=c(iExp,:);
text(x+dx,y+dy,c1,'FontSize',16,'Color',char(colors(iExp)))
hold on

x0=mean(x);
y0=mean(y);
xa=x-x0;
ya=y-y0;
bsyn(iExp)=sum(xa.*ya)/sum(xa.*xa);
U30(iExp)=(-y0+x0*bsyn(iExp))./bsyn(iExp);
r_bsyn(iExp)=min(min(corrcoef(x,y)));
hold on
xt=[min(x)-0.1*(max(x)-min(x)):0.01*(max(x)-min(x)):max(x)+0.1*(max(x)-min(x))];
yt=bsyn(iExp).*(xt-U30(iExp));
plot(xt,yt,char(strcat(colors(iExp),'--')))
disp(char(strcat(Exps(iExp),': bsyn = ',sprintf('%3.3g',bsyn(iExp)),', U30 = ',sprintf('%3.3g',U30(iExp)),'  ;r = ',sprintf('%3.2f',r_bsyn(iExp)))))
end
title('Mid-lat synoptic eddies')
xlabel(strcat('Tropospheric U ',sprintf('%d',lat_Ueady),'N, ',sprintf('%d',lev_Ueady),'mb'))
%ylabel('stationary eddy heat transport')
ylabel(strcat('$\frac{1}{\delta y} \int_{',sprintf('%d',latEPstaint0),'N}^{',sprintf('%d',latEPstaint1),'N} \left.\it{EP}_{\it{s},z} \right|_{',sprintf('%d',levEPsta),'mb} dy $'),'Interpreter','latex')
set(gca,'FontSize',18)
corr(Utropoint(:),EPwint(:))


figure(9)
clf
for iExp=1:NExps
x=Utropo(iExp,:);
y=(EPlk2int(iExp,:)-EP40k2int(iExp,:)-EPstak2int(iExp,:))./deltaytropo;% pure transient
%y=(EPlk2int(iExp,:)-EP40k2int(iExp,:))./deltaytropo;
% y=(EPlk2int(iExp,:)-EPstak2int(iExp,:))./deltaytropo;
% y=(EPint(iExp,:)-EPstaint(iExp,:)-EP40int(iExp,:))./deltaytropo;
% y=EPstak2int(iExp,:)./deltaytropo;
plot(x,y,char(strcat(colors(iExp),'^')))
hold on
dx=0.00002;
dy=0.0001;
c1=c(iExp,:);
text(x+dx,y+dy,c1,'FontSize',16,'Color',char(colors(iExp)))
hold on

x0=mean(x);
y0=mean(y);
xa=x-x0;
ya=y-y0;
b(iExp)=sum(xa.*ya)/sum(xa.*xa);
U300(iExp)=(-y0+x0*b(iExp))./b(iExp);
r_b(iExp)=min(min(corrcoef(x,y)));
hold on
xt=[min(x)-0.1*(max(x)-min(x)):0.01*(max(x)-min(x)):max(x)+0.1*(max(x)-min(x))];
yt=b(iExp).*(xt-U300(iExp));
plot(xt,yt,char(strcat(colors(iExp),'--')))
disp(char(strcat(Exps(iExp),': b = ',sprintf('%3.3g',b(iExp)),', U300 = ',sprintf('%3.3g',U300(iExp)),'  ;r = ',sprintf('%3.2f',r_bsyn(iExp)))))
end
title('Mid-lat large scale transient eddies')
xlabel(strcat('Tropospheric U ',sprintf('%d',lat_Ueady),'N, ',sprintf('%d',lev_Ueady),'mb'))
%ylabel('stationary eddy heat transport')
ylabel(strcat('$\frac{1}{\delta y} \int_{',sprintf('%d',latEPstaint0),'N}^{',sprintf('%d',latEPstaint1),'N} \left.\it{EP}_{l,z} \right|_{',sprintf('%d',levEPsta),'mb} dy $'),'Interpreter','latex')
set(gca,'FontSize',18)

disp(' ')
disp('The correlation between mean jet speed and the mid-lat generated synoptic waves are collapse together.')
disp('========================================================================================')
disp(' ')
pause

%% How well can synoptic eddies represent all other waves excluding MJO-waves
disp('=======================HOW WELL SYNOPTIC EDDIES CAN REPRESENT MID-LAT LONG WAVE=======================')
otherwaves=(EPlk2int-EP40int-EPstak2int)./deltaytropo;
synoptic=EPwint./deltaysta;
figure(21)
clf
for iExp=1:NExps
x=synoptic(iExp,:);
y=otherwaves(iExp,:);
plot(x,y,char(strcat(colors(iExp),'^')))
hold on
dx=0.00002;
dy=0.0001;
c1=c(iExp,:);
text(x+dx,y+dy,c1,'FontSize',16,'Color',char(colors(iExp)))
hold on

x0=mean(x);
y0=mean(y);
xa=x-x0;
ya=y-y0;
alphal_syn(iExp)=sum(xa.*ya)/sum(xa.*xa);
intcpt(iExp)=(-y0+x0*alphal_syn(iExp))./alphal_syn(iExp);
r_alphal(iExp)=min(min(corrcoef(x,y)));
hold on
xt=[min(x)-0.1*(max(x)-min(x)):0.01*(max(x)-min(x)):max(x)+0.1*(max(x)-min(x))];
yt=alphal_syn(iExp).*(xt-intcpt(iExp));
plot(xt,yt,char(strcat(colors(iExp),'--')))
disp(char(strcat(Exps(iExp),': alphal(synoptic) = ',sprintf('%3.3g',alphal_syn(iExp)),', intcpt = ',sprintf('%3.3g',intcpt(iExp)),'  ;r = ',sprintf('%3.2f',r_bsyn(iExp)))))
end
title('Mid-lat generated large-scale eddy ~ synoptic eddy')
xlabel(strcat('synoptic eddies',sprintf('%d',latEPstaint0),'-',sprintf('%d',latEPstaint1),'N'))
%ylabel('stationary eddy heat transport')
ylabel(strcat('large-scale eddies'))
set(gca,'FontSize',18)

disp('Better correlation when boxlat0 is norther (55N).')
disp('===================================================================================')
disp(' ')

%% mid-lat spectra (alphal & alphas)

alphal1=EPlk1int(1,1)/EPstaint(1,1);
alphas1=(EPint(1,1)-EPlk1int(1,1))/EPstaint(1,1);
alphal2=EPlk2int(1,1)/EPstaint(1,1);
alphas2=(EPint(1,1)-EPlk2int(1,1))/EPstaint(1,1);
alphasyn=EPwint(1,1)/EPstaint(1,1);
disp('=======================(BASED ON CONTROL) MID-LAT GENERATED WAVE SPECTRA RATIO=======================')
disp(strcat('alpha_syn=',sprintf('%3.2f',alphasyn)));
disp(strcat('if limit to k=1: alphal=',sprintf('%3.2f',alphal1),', alphas=',sprintf('%3.2f',alphas1)));
disp(strcat('if limit to k=2: alphal=',sprintf('%3.2f',alphal2),', alphas=',sprintf('%3.2f',alphas2)));
disp('===================================================================================')
disp(' ')
pause

%% How mid-lat wind be weakened by absorbing MJO-waves
disp('=======================JET SPEED ~ MJO EP CONVERGENCE=======================')
divEP=(EP40int_T-EP40int_B+EP40int_N-EP40int_S)./(deltaytropo*deltap).*86400;
%divEP=(EPint_T-EPint_B+EPint_N-EPint_S)./(deltaytropo*deltap).*86400;
figure(20)
clf
for iExp=1:NExps
x=divEP(iExp,:);
y=Utropoint(iExp,:);

plot(x,y,char(strcat(colors(iExp),'^')))
hold on
c1=c(iExp,:);
text(x,y,c1,'FontSize',16,'Color',char(colors(iExp)))
hold on

x0=mean(x);
y0=mean(y);
xa=x-x0;
ya=y-y0;
Cu(iExp)=sum(xa.*ya)/sum(xa.*xa);
U20(iExp)=(y0-x0*Cu(iExp));
r_Cu(iExp)=min(min(corrcoef(x,y)));
hold on
xt=[min(x)-0.1*(max(x)-min(x)):0.01*(max(x)-min(x)):max(x)+0.1*(max(x)-min(x))];
yt=Cu(iExp).*xt+U20(iExp);
plot(xt,yt,char(strcat(colors(iExp),'--')))
disp(char(strcat(Exps(iExp),': Cu = ',sprintf('%3.3g',Cu(iExp)),', U20 = ',sprintf('%3.3g',U20(iExp)),'  ;r = ',sprintf('%3.2f',r_Cu(iExp)))))
end
title('Weakened Jet speed by MJO-waves')
ylabel(strcat('Tropospheric U in the box'))
xlabel(strcat('MJO-freq EP convergence in the box'))
set(gca,'FontSize',18)
disp('Better correlation when boxlat0 is souther (40N) and boxlat1 is souther (55N) and boxlev0 is higher (500mb).')
disp('===================================================================================')
disp(' ')
pause

%% reviewer's suggestion: U_xavg
% cs=jet(length(Amps));
% for iExp=1:NExps
%     figure(100+iExp)
%     clf
%     hold on
%     plot(-100,-100,'LineWidth',2,'Color','k','LineStyle','--','Marker','none');
%     plot(-100,-100,'LineWidth',2,'Color','k','LineStyle','-','Marker','none');
%     legend('    U','\Delta U')
%     for icase=1:N
%         %yyaxis left
%             set(gca,'YColor','k')
%             plot(lat,squeeze(U_xavg(iExp,icase,:)-U_xavg(iExp,1,:)),'LineWidth',2,'Color',cs(icase,:),'LineStyle','-','Marker','none')
%             ylim([-20,50])
%             %ylabel(strcat('dU @',sprintf('%d',lat_conv),'N'))
%             %ylabel('\Delta max jet speed (m/s)')
%         %yyaxis right
%             %set(gca,'YColor','k')
%             hold on
%             plot(lat,smooth(squeeze(U_xavg(iExp,icase,:)),10,'lowess'),'LineWidth',2,'Color',cs(icase,:),'LineStyle','--','Marker','none')
%             %ylim([-20,50])
%             %ylabel(strcat('U @',sprintf('%d',lat_conv),'N'))
%             ylabel('Zonally averaged U (m/s)')
%             
%     end
%     xlabel('latitude (^\circE)')
%     xticks([-90:30:90])
%     xlim([-90,90]);
%     title(strcat(char(Exps(iExp)),'Amps'),'Interpreter','none')
%     set(gca,'FontSize',18)
%     set(gcf,'Position',[684   517   601   369])
%    
% end
end