clc
clear
close all;

load('F:/OHC_trend_by_IAP/CESM_xhgh/OMEGA/CESM_xghg_Q_1000.mat')
load('F:/OHC_trend_by_IAP/CESM_xhgh/OMEGA/CESM_xghg_OMEGA_500.mat')

load('F:/OHC_trend_by_IAP/CESM_xaer/OMEGA/CESM_xaer_Q_1000.mat')
load('F:/OHC_trend_by_IAP/CESM_xaer/OMEGA/CESM_xaer_OMEGA_500.mat')
%--------------------------------------------------------------------------
%---------xghg------------------------------
Q_xghg_climate_1000=squeeze(nanmean(Q_xghg_1000_yr_r(:,:,:,:),3));
OMEGA_xghg_climate_500=squeeze(nanmean(OMEGA_xghg_500_yr_r(:,:,:,:),3));


Q_xghg_a=squeeze(nanmean(Q_xghg_1000_yr_r(:,:,(1985-1920)+1:(2005-1920)+1,:),3))-...
    squeeze(nanmean(Q_xghg_1000_yr_r(:,:,(1920-1920)+1:(1940-1920)+1,:),3));

OMEGA_xghg_a=squeeze(nanmean(OMEGA_xghg_500_yr_r(:,:,(1985-1920)+1:(2005-1920)+1,:),3))-...
    squeeze(nanmean(OMEGA_xghg_500_yr_r(:,:,(1920-1920)+1:(1940-1920)+1,:),3));

%Dynamic 
PRE_xghg_do=Q_xghg_climate_1000.*OMEGA_xghg_a;
%Themodynamic
PRE_xghg_re=OMEGA_xghg_climate_500.*Q_xghg_a;

%model mean
PRE_xghg_do_modelmean=squeeze(nanmean(PRE_xghg_do,3));
PRE_xghg_re_modelmean=squeeze(nanmean(PRE_xghg_re,3));
% %T-test
% [nxx nyy nrr]=size(PRE_xghg_do);
% for i=1:nxx
%     for j=1:nyy
%         PRE_xghg_do_tt(i,j)=ttest(PRE_xghg_do(i,j,:),0,'Alpha',0.05);
%          PRE_xghg_re_tt(i,j)=ttest(PRE_xghg_re(i,j,:),0,'Alpha',0.05);
%     end
% end
% 
% PRE_xghg_do_modelmean_new=PRE_xghg_do_modelmean;
% PRE_xghg_do_modelmean_new(PRE_xghg_do_tt==0)=NaN;
% PRE_xghg_re_modelmean_new=PRE_xghg_re_modelmean;
% PRE_xghg_re_modelmean_new(PRE_xghg_re_tt==0)=NaN;
%-----------xaer-----------------------------------------------------------
Q_xaer_climate_1000=squeeze(nanmean(Q_xaer_1000_yr_r(:,:,:,:),3));
OMEGA_xaer_climate_500=squeeze(nanmean(OMEGA_xaer_500_yr_r(:,:,:,:),3));

Q_xaer_a=squeeze(nanmean(Q_xaer_1000_yr_r(:,:,(1985-1920)+1:(2005-1920)+1,:),3))-...
    squeeze(nanmean(Q_xaer_1000_yr_r(:,:,(1920-1920)+1:(1940-1920)+1,:),3));

OMEGA_xaer_a=squeeze(nanmean(OMEGA_xaer_500_yr_r(:,:,(1985-1920)+1:(2005-1920)+1,:),3))-...
    squeeze(nanmean(OMEGA_xaer_500_yr_r(:,:,(1920-1920)+1:(1940-1920)+1,:),3));

%Dynamic 
PRE_xaer_do=Q_xaer_climate_1000.*OMEGA_xaer_a;
%Themodynamic
PRE_xaer_re=OMEGA_xaer_climate_500.*Q_xaer_a;

%model mean
PRE_xaer_do_modelmean=squeeze(nanmean(PRE_xaer_do,3));
PRE_xaer_re_modelmean=squeeze(nanmean(PRE_xaer_re,3));
% %T-test
% [nxx nyy nrr]=size(PRE_xghg_do);
% for i=1:nxx
%     for j=1:nyy
%         PRE_xaer_do_tt(i,j)=ttest(PRE_xaer_do(i,j,:),0,'Alpha',0.05);
%          PRE_xaer_re_tt(i,j)=ttest(PRE_xaer_re(i,j,:),0,'Alpha',0.05);
%     end
% end
% 
% PRE_xaer_do_modelmean_new=PRE_xaer_do_modelmean;
% PRE_xaer_do_modelmean_new(PRE_xaer_do_tt==0)=NaN;
% PRE_xaer_re_modelmean_new=PRE_xaer_re_modelmean;
% PRE_xaer_re_modelmean_new(PRE_xaer_re_tt==0)=NaN;

%----his-------------------------------------------------------------------
load('F:/IO_Dome/CESM_his/Q/CESM_his_Q_1000.mat')
load('F:/IO_Dome/CESM_his/OMEGA/CESM_his_OMEGA_500.mat')

latmin=-40;
latmax=40;
lonmin=0;
lonmax=360;

Q_climate_1000=squeeze(nanmean(Q_1000_yr_r(:,:,:,1:34),3));
Nlat=find(LAT>=latmin&LAT<=latmax);
OMEGA_climate_500=squeeze(nanmean(OMEGA_500_yr_r(:,Nlat,:,:),3));

Q_a=squeeze(nanmean(Q_1000_yr_r(:,:,(1985-1920)+1:(2005-1920)+1,1:34),3))-...
    squeeze(nanmean(Q_1000_yr_r(:,:,(1920-1920)+1:(1940-1920)+1,1:34),3));

OMEGA_a=squeeze(nanmean(OMEGA_500_yr_r(:,Nlat,(1985-1920)+1:(2005-1920)+1,:),3))-...
    squeeze(nanmean(OMEGA_500_yr_r(:,Nlat,(1920-1920)+1:(1940-1920)+1,:),3));

%Dynamic 
PRE_do=Q_climate_1000.*OMEGA_a;
%Themodynamic
PRE_re=OMEGA_climate_500.*Q_a;

%model mean
PRE_do_modelmean=squeeze(nanmean(PRE_do,3));
PRE_re_modelmean=squeeze(nanmean(PRE_re,3));

% %T-test
% [nxx nyy nrr]=size(PRE_do);
% for i=1:nxx
%     for j=1:nyy
%         PRE_do_tt(i,j)=ttest(PRE_do(i,j,:),0,'Alpha',0.05);
%          PRE_re_tt(i,j)=ttest(PRE_re(i,j,:),0,'Alpha',0.05);
%     end
% end
% 
% PRE_do_modelmean_new=PRE_do_modelmean;
% PRE_do_modelmean_new(PRE_do_tt==0)=NaN;
% PRE_re_modelmean_new=PRE_re_modelmean;
% PRE_re_modelmean_new(PRE_re_tt==0)=NaN;
%% -------------------------------------------------------------------------
PRE_aer_do_modelmean=PRE_do_modelmean-PRE_xaer_do_modelmean;
PRE_aer_re_modelmean=PRE_re_modelmean-PRE_xaer_re_modelmean;

PRE_ghg_do_modelmean=PRE_do_modelmean-PRE_xghg_do_modelmean;
PRE_ghg_re_modelmean=PRE_re_modelmean-PRE_xghg_re_modelmean;

[X,Y]=meshgrid(LON_q,LAT_q);
%% ----------------TS trend------------
load('H:\SCS_SSS\图9 降水分析single\CESM_his_TS.mat')
latmin=-40;
latmax=40;
lonmin=0;
lonmax=360;

Nlat=find(Lat>=latmin&Lat<=latmax);
TS_his_a_1=TS_his_a(:,Nlat);
TS_aer_a_1=TS_aer_a(:,Nlat);
TS_ghg_a_1=TS_ghg_a(:,Nlat);
%% %-----------------fugure---------------------------
figure('Units','centimeters','Position',[8, 2, 21, 18]);
subplot('position',[0.55 0.7 0.42 0.25]);
contourf(X,Y,-PRE_do_modelmean'*10^5,30,'linestyle','none');
hold on
[c h]=contour(X,Y,TS_his_a_1',[0: 0.25:2],'-k','LineWidth',1.5);
clabel(c,h,'color',[0 0 0],'fontsize',10);
[c h]=contour(X,Y,TS_his_a_1',[-2: 0.25:0],'--k','LineWidth',1);
clabel(c,h,'color',[0 0 0],'fontsize',10);
grid on
load coast;
long=long+(long<0)*360;
long(long>359)=NaN;
plot(long,lat,'color','k')
ylim([-40 40])
xlim([30 300])
caxis([-15 15]);
set(gca,'xtick',30:60:300)
set(gca,'Xticklabel',{'30E','90E','150E','180',...
    '150W','90W','30W'})
set(gca,'Ytick',-60:20:60,'Yticklabel',{'60S','40S','20S','EQ',...
    '20N','40N','60N'},'fontname','Times New Roman','fontsize',10,'fontweight','normal')
title('(b) CESM-his Dynamic Contributions','fontsize',12,'FontName','Times New Roman')

load('F:/SCS_SSS/图6降水分析/CESM_his_PRECL_climate.mat')
subplot('position',[0.05 0.7 0.42 0.25]);
contourf(X,Y,-PRE_re_modelmean'*10^5,30,'linestyle','none');
hold on
[c h]=contour(X,Y,PRE_climate','color',[0 0 0],'linewidth',1.5);
clabel(c,h,'color',[0 0 0],'fontsize',10);
grid on
load coast;
long=long+(long<0)*360;
long(long>359)=NaN;
plot(long,lat,'color','k')
ylim([-40 40])
xlim([30 300])
caxis([-3 3]);
set(gca,'xtick',30:60:300)
set(gca,'Xticklabel',{'30E','90E','150E','180',...
    '150W','90W','30W'})
set(gca,'Ytick',-60:20:60,'Yticklabel',{'60S','40S','20S','EQ',...
    '20N','40N','60N'},'fontname','Times New Roman','fontsize',10,'fontweight','normal')
title('(a) CESM-his Thermodynamic Contributions','fontsize',12,'FontName','Times New Roman')

subplot('position',[0.55 0.39 0.42 0.25]);
contourf(X,Y,-PRE_aer_do_modelmean'*10^5,30,'linestyle','none');
hold on
[c h]=contour(X,Y,TS_aer_a_1',[0: 0.25:2],'-k','LineWidth',1.5);
clabel(c,h,'color',[0 0 0],'fontsize',10);
[c h]=contour(X,Y,TS_aer_a_1',[-2: 0.25:0],'--k','LineWidth',1);
clabel(c,h,'color',[0 0 0],'fontsize',10);
grid on
load coast;
long=long+(long<0)*360;
long(long>359)=NaN;
plot(long,lat,'color','k')
ylim([-40 40])
xlim([30 300])
caxis([-15 15]);
set(gca,'xtick',30:60:300)
set(gca,'Xticklabel',{'30E','90E','150E','180',...
    '150W','90W','30W'})
set(gca,'Ytick',-60:20:60,'Yticklabel',{'60S','40S','20S','EQ',...
    '20N','40N','60N'},'fontname','Times New Roman','fontsize',10,'fontweight','normal')
title('(d) CESM-Aer Dynamic Contributions','fontsize',12,'FontName','Times New Roman')


subplot('position',[0.05 0.39 0.42 0.25]);
contourf(X,Y,-PRE_aer_re_modelmean'*10^5,30,'linestyle','none');
hold on
[c h]=contour(X,Y,PRE_climate','color',[0 0 0],'linewidth',1.5);
clabel(c,h,'color',[0 0 0],'fontsize',10);
grid on
load coast;
long=long+(long<0)*360;
long(long>359)=NaN;
plot(long,lat,'color','k')
ylim([-40 40])
xlim([30 300])
caxis([-3 3]);
set(gca,'xtick',30:60:300)
set(gca,'Xticklabel',{'30E','90E','150E','180',...
    '150W','90W','30W'})
set(gca,'Ytick',-60:20:60,'Yticklabel',{'60S','40S','20S','EQ',...
    '20N','40N','60N'},'fontname','Times New Roman','fontsize',10,'fontweight','normal')
% m_grid('xtick',6,'ytick',4,'tickdir','out','yaxislocation','left','fontsize',10);%grid璁剧疆
% m_coast('color','k','linewidth',1.5);
title('(c) CESM-Aer Thermodynamic Contributions','fontsize',12,'FontName','Times New Roman')
%----------------ghg------------------------------------------
subplot('position',[0.55 0.08 0.42 0.25]);
contourf(X,Y,-PRE_ghg_do_modelmean'*10^5,30,'linestyle','none');
hold on
grid on
[c h]=contour(X,Y,TS_ghg_a_1',[0: 0.25:2],'-k','LineWidth',1.5);
clabel(c,h,'color',[0 0 0],'fontsize',10);
[c h]=contour(X,Y,TS_ghg_a_1',[-2: 0.25:0],'--k','LineWidth',1.5);
clabel(c,h,'color',[0 0 0],'fontsize',10);
load coast;
long=long+(long<0)*360;
long(long>359)=NaN;
plot(long,lat,'color','k')
ylim([-40 40])
xlim([30 300])
caxis([-15 15]);
set(gca,'xtick',30:60:300)
set(gca,'Xticklabel',{'30E','90E','150E','180',...
    '150W','90W','30W'})
set(gca,'Ytick',-60:20:60,'Yticklabel',{'60S','40S','20S','EQ',...
    '20N','40N','60N'},'fontname','Times New Roman','fontsize',10,'fontweight','normal')
title('(f) CESM-GHGs Dynamic Contributions','fontsize',12,'FontName','Times New Roman')

load('H:/color/Green_yellow.mat')
L=length(colors);
L2=floor(L/2);
a=colors(1:L2,:);b=colors(L2:end-10,:);
c=[a;b];
colormap(c);
hcb=colorbar('location','SouthOutside','position',[0.55,0.03,0.4,0.016]);

subplot('position',[0.05 0.08 0.42 0.25]);
contourf(X,Y,-PRE_ghg_re_modelmean'*10^5,30,'linestyle','none');
hold on
[c h]=contour(X,Y,PRE_climate','color',[0 0 0],'linewidth',1.5);
clabel(c,h,'color',[0 0 0],'fontsize',10);
hold on
grid on
load coast;
long=long+(long<0)*360;
long(long>359)=NaN;
plot(long,lat,'color','k')
ylim([-40 40])
xlim([30 300])
caxis([-3 3]);
set(gca,'xtick',30:60:300)
set(gca,'Xticklabel',{'30E','90E','150E','180',...
    '150W','90W','30W'})
set(gca,'Ytick',-60:20:60,'Yticklabel',{'60S','40S','20S','EQ',...
    '20N','40N','60N'},'fontname','Times New Roman','fontsize',10,'fontweight','normal')
title('(e) CESM-GHGs Thermodynamic Contributions','fontsize',12,'FontName','Times New Roman')
hcb=colorbar('location','SouthOutside','position',[0.05,0.03,0.4,0.016]);
% print(gcf,'-dpng','-r600','CESM_Q_100_single(ghg)');
% print(gcf,'-depsc2','-r600','CESM_Q_100_single(ghg)')
% print -depsc2 -painters CESM_Q_100_single.eps