clc
clear
close all;

load('I:/OHC_trend_by_IAP/CESM_xaer/OMEGA/CESM_xaer_OMEGA_500.mat')
load('I:/OHC_trend_by_IAP/CESM_xbmb/OMEGA/CESM_xbmb_OMEGA_500.mat')
load('I:/OHC_trend_by_IAP/CESM_xhgh/OMEGA/CESM_xghg_OMEGA_500.mat')

load('I:/IO_Dome/CESM_his/Q/CESM_his_Q_1000.mat')
load('I:/IO_Dome/CESM_his/OMEGA/CESM_his_OMEGA_500.mat')

latmin=-40;
latmax=40;
Nlat=find(LAT>=latmin&LAT<=latmax);
OMEGA_his_500_yr_r=squeeze(OMEGA_500_yr_r(:,Nlat,:,:));



OMEGA_aer_500=squeeze(nanmean(OMEGA_his_500_yr_r,4))-squeeze(nanmean(OMEGA_xaer_500_yr_r,4));
OMEGA_ghg_500=squeeze(nanmean(OMEGA_his_500_yr_r,4))-squeeze(nanmean(OMEGA_xghg_500_yr_r,4));
OMEGA_bmb_500=squeeze(nanmean(OMEGA_his_500_yr_r,4))-squeeze(nanmean(OMEGA_xbmb_500_yr_r,4));

OMEGA_aer_500_1=OMEGA_aer_500(:,:,(1920-1920)+1:(1980-1920)+1);
OMEGA_ghg_500_1=OMEGA_ghg_500(:,:,(1920-1920)+1:(1980-1920)+1);
OMEGA_his_500_1=squeeze(nanmean(OMEGA_his_500_yr_r(:,:,(1920-1920)+1:(1980-1920)+1),4));
%trend
[nxx nyy ntt]=size(OMEGA_aer_500_1);
t2=[1:1:ntt];
for i=1:nxx
    for j=1:nyy
        y=squeeze(OMEGA_his_500_1(i,j,:));
        P=polyfit(t2,y',1);
        OMEGA_his_trend(i,j)=P(1);
        
        y=squeeze(OMEGA_aer_500_1(i,j,:));
        P=polyfit(t2,y',1);
        OMEGA_aer_500_trend(i,j)=P(1);

        y=squeeze(OMEGA_ghg_500_1(i,j,:));
        P=polyfit(t2,y',1);
        OMEGA_ghg_500_trend(i,j)=P(1);

%         y=squeeze(OMEGA_bmb_500(i,j,:));
%         P=polyfit(t2,y',1);
%         OMEGA_bmb_500_trend(i,j)=P(1);
    end
end


%% wind
load('I:/OHC_trend_by_IAP/CESM_xaer/wind/CESM_xaer_Vonp_850.mat');
load('I:/OHC_trend_by_IAP/CESM_xhgh/wind/CESM_xghg_Vonp_850.mat')
load('CESM_his_Uonp_850_trend.mat')

VonP_aer_850_trend = VonP_850_trend-VonP_xaer_850_trend;
UonP_aer_850_trend = UonP_850_trend-UonP_xaer_850_trend;
VonP_ghg_850_trend = VonP_850_trend-VonP_xghg_850_trend;
UonP_ghg_850_trend = UonP_850_trend-UonP_xghg_850_trend;

% %MK检验
% for xi=1:size(OMEGA_ghg_500_1,1)
%     for yi=1:size(OMEGA_ghg_500_1,2)
%         OMEGA_ghg=squeeze(OMEGA_ghg_500_1(xi,yi,:));
%         [H(xi,yi),p_value]=Mann_Kendall(OMEGA_ghg,0.1);
% 
%         OMEGA_are=squeeze(OMEGA_aer_500_1(xi,yi,:));
%         [H2(xi,yi),p_value]=Mann_Kendall(OMEGA_are,0.1);
% 
%         OMEGA_his=squeeze(OMEGA_his_500_1(xi,yi,:));
%         [H3(xi,yi),p_value]=Mann_Kendall(OMEGA_his,0.1);
%     end
% end
% H=double(H);H2=double(H2);H3=double(H3);
VonP_850_trend(abs(VonP_850_trend)<=10^-3)=NaN;UonP_850_trend(abs(UonP_850_trend)<=10^-3)=NaN;
VonP_ghg_850_trend(abs(VonP_ghg_850_trend)<=10^-3)=NaN;UonP_ghg_850_trend(abs(UonP_ghg_850_trend)<=10^-3)=NaN;
VonP_aer_850_trend(abs(VonP_aer_850_trend)<=10^-3)=NaN;UonP_aer_850_trend(abs(UonP_aer_850_trend)<=10^-3)=NaN;
%%
figure('Units','centimeters','Position',[8,2,21,15]);
subplot('position',[0.07 0.71 0.8 0.25]);
[long,latt]=meshgrid(double(LON_q),double(LAT_q));
contourf(long,latt,-OMEGA_his_trend'*12*100,500,'linestyle','none');
hold on;
% quiver(long(1:4:end,1:4:end),latt(1:4:end,1:4:end),...
%     UonP_850_trend(1:4:end,1:4:end)'*1200,VonP_850_trend(1:4:end,1:4:end)'*1200,0,'Color', 'k');
load coast;
long=long+(long<0)*360;
long(long>359)=NaN;
plot(long,lat,'color','k')
ylim([-40 40])
xlim([0 360])
set(gca,'xtick',0:30:360)
set(gca,'Xticklabel',{'0','30E','60E','90E','120E','150E','180',...
    '150W','120W','90W','60W','30W'})
set(gca,'Ytick',-40:20:40,'Yticklabel',{'40S','20S','EQ',...
    '20N','40N'},'fontname','Times New Roman','fontsize',10,'fontweight','normal')
% m_grid('xtick',12,'ytick',7,'tickdir','out','yaxislocation','left','fontsize',10);%grid璁剧疆
caxis([-0.21 0.21]);
title('(a) CESM-his OMEGA 500hPa & wind 850hPa trend','fontsize',12)


subplot('position',[0.07 0.05 0.8 0.25]);
[long,latt]=meshgrid(double(LON_q),double(LAT_q));
contourf(long,latt,-OMEGA_ghg_500_trend'*12*100,500,'linestyle','none');
hold on;
% quiver(long(1:4:end,1:4:end),latt(1:4:end,1:4:end),...
%     UonP_ghg_850_trend(1:4:end,1:4:end)'*1200,VonP_ghg_850_trend(1:4:end,1:4:end)'*1200,0,'Color', 'k');
load coast;
long=long+(long<0)*360;
long(long>359)=NaN;
plot(long,lat,'color','k')
ylim([-40 40])
xlim([0 360])
set(gca,'xtick',0:30:360)
set(gca,'Xticklabel',{'0','30E','60E','90E','120E','150E','180',...
    '150W','120W','90W','60W','30W'})
set(gca,'Ytick',-40:20:40,'Yticklabel',{'40S','20S','EQ',...
    '20N','40N'},'fontname','Times New Roman','fontsize',10,'fontweight','normal')
% m_grid('xtick',12,'ytick',7,'tickdir','out','yaxislocation','left','fontsize',10);%grid璁剧疆
caxis([-0.21 0.21]);
title('(c) CESM-GHGs OMEGA 500hPa & wind 850hPa trend','fontsize',12)

% subplot(2,1,2)
subplot('position',[0.07 0.38 0.8 0.25]);
[long,latt]=meshgrid(double(LON_q),double(LAT_q));
contourf(long,latt,-OMEGA_aer_500_trend'*12*100,500,'linestyle','none');
hold on;
% quiver(long(1:4:end,1:4:end),latt(1:4:end,1:4:end),...
%     UonP_aer_850_trend(1:4:end,1:4:end)'*1200,VonP_aer_850_trend(1:4:end,1:4:end)'*1200,0,'Color', 'k');
quiver(30,30,...
    10',0',0,'Color', 'r');
load coast;
long=long+(long<0)*360;
long(long>359)=NaN;
plot(long,lat,'color','k')
ylim([-40 40])
xlim([0 360])
set(gca,'xtick',0:30:360)
set(gca,'Xticklabel',{'0','30E','60E','90E','120E','150E','180',...
    '150W','120W','90W','60W','30W'})
set(gca,'Ytick',-40:20:40,'Yticklabel',{'40S','20S','EQ',...
    '20N','40N'},'fontname','Times New Roman','fontsize',10,'fontweight','normal')
% m_grid('xtick',12,'ytick',7,'tickdir','out','yaxislocation','left','fontsize',10);%grid璁剧疆
caxis([-0.21 0.21]);
title('(b) CESM-Aer OMEGA 500hPa & wind 850hPa trend','fontsize',12)

load('H:/color/b_r3.mat')
L=length(colors);
L2=floor(L/2);
a=colors(1:L2,:);b=colors(L2:end,:);
c=[a;b];
colormap(c);
hcb=colorbar('position',[0.92,0.17,0.015,0.75]);
% print(gcf,'-dpng','-r600','CESM_OMEGA_single_trend(aer500hpa_1920_1980)');
% print(gcf,'-depsc2','-r600','CESM_OMEGA_single_trend(aer500hpa)')
% print -depsc2 -painters CESM_OMEGA_single_trend(aer500hpa).eps
%%

% %trend
% [nxx nyy ntt nrr]=size(OMEGA_his_500_yr_r);
% t2=[1:1:ntt];
% for r=1:nrr
%     for i=1:nxx
%         for j=1:nyy
%             y=squeeze(OMEGA_his_500_yr_r(i,j,:,r));
%             P=polyfit(t2,y',1);
%             OMEGA_his_500_trend(i,j,r)=P(1);
%         end
%     end
% end
% [nxx nyy nrr]=size(OMEGA_his_500_trend);
% for i=1:nxx
%     for j=1:nyy
%         OMEGA_his_500_trend_tt(i,j)=ttest(OMEGA_his_500_trend(i,j,:),0,'Alpha',0.05);
%     end
% end
% OMEGA_his_trend=squeeze(nanmean(OMEGA_his_500_trend,3));
% 
% OMEGA_his_trend_new=OMEGA_his_trend;
% OMEGA_his_trend_new(OMEGA_his_500_trend_tt==0)=NaN;
% %trend
% [nxx nyy ntt nrr]=size(OMEGA_xaer_500_yr_r);
% t2=[1:1:ntt];
% for r=1:nrr
%     for i=1:nxx
%         for j=1:nyy
%             y=squeeze(OMEGA_xaer_500_yr_r(i,j,:,r));
%             P=polyfit(t2,y',1);
%             OMEGA_xaer_500_trend(i,j,r)=P(1);
% 
%             y=squeeze(OMEGA_xghg_500_yr_r(i,j,:,r));
%             P=polyfit(t2,y',1);
%             OMEGA_xghg_500_trend(i,j,r)=P(1);
%         end
%     end
% end
% [nxx nyy nrr]=size(OMEGA_his_500_trend);
% for i=1:nxx
%     for j=1:nyy
%         OMEGA_his_500_trend_tt(i,j)=ttest(OMEGA_his_500_trend(i,j,:),0,'Alpha',0.05);
%     end
% end
% OMEGA_his_trend=squeeze(nanmean(OMEGA_his_500_trend,3));
% 
% OMEGA_his_trend_new=OMEGA_his_trend;
% OMEGA_his_trend_new(OMEGA_his_500_trend_tt==0)=NaN;
% 
% 
% %modelmean
% OMEGA_his_trend=squeeze(nanmean(OMEGA_his_500_trend,3));
% OMEGA_xaer_trend=squeeze(nanmean(OMEGA_xaer_500_trend,3));
% OMEGA_xghg_trend=squeeze(nanmean(OMEGA_xghg_500_trend,3));
% 
% OMEGA_aer_trend_2=OMEGA_his_trend-OMEGA_xaer_trend;
% OMEGA_ghg_trend_2=OMEGA_his_trend-OMEGA_xghg_trend;

% % figure('Units','centimeters','Position',[8,2,21,15]);
% subplot('position',[0.07 0.71 0.8 0.25]);
% % m_proj('Mercator','lat',[min(LAT_q) max(LAT_q)],'lon',[min(LON_q) max(LON_q)]);
% [long,latt]=meshgrid(double(LON_q),double(LAT_q));
% contourf(long,latt,-OMEGA_his_trend',500,'linestyle','none');
% % m_coast('patch',[.6 .6 .6],'linewidth',1)
% %  m_coast('color','k','linewidth',1.5);
% hold on;
% load coast;
% long=long+(long<0)*360;
% long(long>359)=NaN;
% plot(long,lat,'color','k')
% ylim([-40 40])
% xlim([0 360])
% set(gca,'xtick',0:30:360)
% set(gca,'Xticklabel',{'0','30E','60E','90E','120E','150E','180',...
%     '150W','120W','90W','60W','30W'})
% set(gca,'Ytick',-40:20:40,'Yticklabel',{'40S','20S','EQ',...
%     '20N','40N'},'fontname','Times New Roman','fontsize',10,'fontweight','normal')
% % m_grid('xtick',12,'ytick',7,'tickdir','out','yaxislocation','left','fontsize',10);%grid璁剧疆
% caxis([-13*10^-5 13*10^-5]);
% title('(a) CESM-his OMEGA 500hPa trend','fontsize',12)
% 
% load('/Volumes/wu/color/r_w_b3.mat')
% L=length(colors);
% L2=floor(L/2);
% a=colors(1:L2,:);b=colors(L2:end,:);
% c=[a;b];
% colormap(c);
% hcb=colorbar('location','SouthOutside','position',[0.3,0.15,0.45,0.015]);
% % print(gcf,'-dpng','-r600','CESM_OMEGA_single_trend(his500hpa)');
% % print(gcf,'-depsc2','-r600','CESM_OMEGA_single_trend(his500hpa)')
% % print -depsc2 -painters CESM_OMEGA_single_trend(his500hpa).eps
% print -depsc2 -painters CESM_OMEGA_single_trend(aer500hpa).eps