clc
clear
close all;
load('H:\SCS_SSS\图4aer_ghg_SSS\CESM_S_0_55_single.mat')
load('F:\OHC_trend_by_IAP\CESM_xhgh\CESM_xhgh_SCS_0_55_trend.mat')
load('F:\OHC_trend_by_IAP\CESM_xaer\CESM_xaer_SCS_0_55_trend.mat')
load('F:\OHC_trend_by_IAP\CESM_xbmb\CESM_xbmb_SCS_0_55_trend.mat')

SSS_his=squeeze(nanmean(S_0_55_r(:,:,(1920-1920)+1:(2005-1920)+1,:),4));
SSS_xaer=squeeze(nanmean(S_xaer_0_55_r(:,:,(1920-1920)+1:(2005-1920)+1,:),4));
SSS_xhgh=squeeze(nanmean(S_xghg_0_55_r(:,:,(1920-1920)+1:(2005-1920)+1,:),4));
SSS_xbmb=squeeze(nanmean(S_xbmb_0_55_r(:,:,(1920-1920)+1:(2005-1920)+1,:),4));
%
SSS_aer=SSS_his-SSS_xaer;
SSS_hgh=SSS_his-SSS_xhgh;
SSS_bmb=SSS_his-SSS_xbmb;

% trend
[nyy nzz ntt]=size(SSS_aer);
t2=[1:1:ntt];
for i=1:nyy
    for j=1:nzz
        y=squeeze(SSS_aer(i,j,:));
        P=polyfit(t2,y',1);
        SSS_aer_trend(i,j)=P(1);

        y=squeeze(SSS_hgh(i,j,:));
        P=polyfit(t2,y',1);
        SSS_hgh_trend(i,j)=P(1);

        y=squeeze(SSS_bmb(i,j,:));
        P=polyfit(t2,y',1);
        SSS_bmb_trend(i,j)=P(1);
    end
end

%MK检验
for xi=1:size(SSS_hgh,1)
    for yi=1:size(SSS_hgh,2)
        salt=squeeze(SSS_hgh(xi,yi,:));
        [H(xi,yi),p_value]=Mann_Kendall(salt,0.05);
    end
end
H=double(H);

for xi=1:size(SSS_aer,1)
    for yi=1:size(SSS_aer,2)
        salt=squeeze(SSS_aer(xi,yi,:));
        [H2(xi,yi),p_value]=Mann_Kendall(salt,0.05);
    end
end
H2=double(H2);

for xi=1:size(SSS_bmb,1)
    for yi=1:size(SSS_bmb,2)
        salt=squeeze(SSS_bmb(xi,yi,:));
        [H3(xi,yi),p_value]=Mann_Kendall(salt,0.05);
    end
end
H3=double(H3);

figure('Units','centimeters','Position',[8,2,21,15]);
subplot(2,1,1)
m_proj('Mercator','lat',[min(LAT) max(LAT)],'lon',[min(LON) max(LON)]);
[long,latt]=meshgrid(double(LON),double(LAT));
m_contourf(long,latt,SSS_hgh_trend'*1200,500,'linestyle','none');
m_coast('patch',[.7 .7 .7],'linewidth',1,'edgecolor',[.5 .5 .5])
hold on;
[nxx nyy ntt]=size(SSS_hgh_trend);
  for i=1:nxx
     for j=1:nyy
         if(H(i,j)==1)
             m_plot(long(j,i),latt(j,i),'k.','markersize',3.5);
         end
%          UE(1:2:end,1:2:end)=NaN;
     end
  end  
m_grid('xtick',10,'ytick',7,'tickdir','out','yaxislocation','left','fontsize',10);%grid设置
caxis([-5 5]);
% rectangle('Position',[-0.58,-0.07,0.32,0.5],'LineWidth',2,'LineStyle','-');
% lon_x=109:1:120; lat_y=0.8*lon_x-85.5;
% m_plot(lon_x,lat_y,'k','linewidth',2);
lon_x1=101:120;lat_y1=ones(size(lon_x1))*1.7;
m_plot(lon_x1,lat_y1,'k','linewidth',2)
lat_y2=1.7:0.1:25.5;lon_x2=ones(size(lat_y2))*101;
m_plot(lon_x2,lat_y2,'k','linewidth',2)
lat_y3=1.7:0.1:25.5;lon_x3=ones(size(lat_y3))*120;
m_plot(lon_x3,lat_y3,'k','linewidth',2)
lon_x5=101:120;lat_y5=ones(size(lon_x5))*25.5;
m_plot(lon_x5,lat_y5,'k','linewidth',2)
title('(a)CESM-GHG S0-55m trend','fontsize',15)

% figure('Units','centimeters','Position',[8,2,21,15]);
subplot(2,1,2)
m_proj('Mercator','lat',[min(LAT) max(LAT)],'lon',[min(LON) max(LON)]);
[long,latt]=meshgrid(double(LON),double(LAT));
m_contourf(long,latt,SSS_aer_trend'*1200,500,'linestyle','none');
m_coast('patch',[.7 .7 .7],'linewidth',1,'edgecolor',[.5 .5 .5])
hold on;
[nxx nyy ntt]=size(SSS_hgh_trend);
  for i=1:nxx
     for j=1:nyy
         if(H2(i,j)==1)
             m_plot(long(j,i),latt(j,i),'k.','markersize',3.5);
         end
%          UE(1:2:end,1:2:end)=NaN;
     end
  end 
m_grid('xtick',10,'ytick',7,'tickdir','out','yaxislocation','left','fontsize',10);%grid设置
caxis([-5 5]);
% rectangle('Position',[-0.58,-0.07,0.32,0.5],'LineWidth',2,'LineStyle','-');
% lon_x=109:1:120; lat_y=0.8*lon_x-85.5;
% m_plot(lon_x,lat_y,'k','linewidth',2);
lon_x1=101:120;lat_y1=ones(size(lon_x1))*1.7;
m_plot(lon_x1,lat_y1,'k','linewidth',2)
lat_y2=1.7:0.1:25.5;lon_x2=ones(size(lat_y2))*101;
m_plot(lon_x2,lat_y2,'k','linewidth',2)
lat_y3=1.7:0.1:25.5;lon_x3=ones(size(lat_y3))*120;
m_plot(lon_x3,lat_y3,'k','linewidth',2)
lon_x5=101:120;lat_y5=ones(size(lon_x5))*25.5;
m_plot(lon_x5,lat_y5,'k','linewidth',2)
title('(c)CESM-Aer S0-55m trend','fontsize',15)

load('J:/盐度/年代际变化/color/b_r3')
L=length(colors);
L2=floor(L/2);
a=colors(1:L2-5,:);b=colors(L2+5:end,:);
c=[a;b];
colormap(c);
hcb=colorbar('position',[0.85,0.15,0.015,0.75]);
% print(gcf,'-dpng','-r600','CESM_S_0_55_single');
% print(gcf,'-depsc2','-r600','CESM_S_0_55_single')
% print -depsc2 -painters CESM_S_0_55_single.eps
%% ----------------------------------------------------------------------
figure('Units','centimeters','Position',[8,2,21,15]);
m_proj('Mercator','lat',[min(LAT) max(LAT)],'lon',[min(LON) max(LON)]);
[long,latt]=meshgrid(double(LON),double(LAT));
m_contourf(long,latt,SSS_bmb_trend',500,'linestyle','none');
m_coast('patch',[.6 .6 .6],'linewidth',1)
hold on;
[nxx nyy ntt]=size(SSS_bmb_trend);
  for i=1:nxx
     for j=1:nyy
         if(H3(i,j)==1)
             m_plot(long(j,i),latt(j,i),'k.','markersize',3.5);
         end
%          UE(1:2:end,1:2:end)=NaN;
     end
  end 
m_grid('xtick',10,'ytick',7,'tickdir','out','yaxislocation','left','fontsize',10);%grid设置
caxis([-8*10^-3 8*10^-3]);
% rectangle('Position',[-0.58,-0.07,0.32,0.5],'LineWidth',2,'LineStyle','-');
% lon_x=109:1:120; lat_y=0.8*lon_x-85.5;
% m_plot(lon_x,lat_y,'k','linewidth',2);
lon_x1=101:120;lat_y1=ones(size(lon_x1))*1.7;
m_plot(lon_x1,lat_y1,'k','linewidth',2)
lat_y2=1.7:0.1:25.5;lon_x2=ones(size(lat_y2))*101;
m_plot(lon_x2,lat_y2,'k','linewidth',2)
lat_y3=1.7:0.1:25.5;lon_x3=ones(size(lat_y3))*120;
m_plot(lon_x3,lat_y3,'k','linewidth',2)
lon_x5=101:120;lat_y5=ones(size(lon_x5))*25.5;
m_plot(lon_x5,lat_y5,'k','linewidth',2)
title('CESM-bmb S0-55m 1920-2005 trend','fontsize',15)

load('J:/盐度/年代际变化/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.15,0.015,0.75]);
% print(gcf,'-dpng','-r600','CESM_S_0_55_single(bmb)');
%print -depsc2 -painters CESM_S_0_55_single(time_line)_bmb.eps
%%
%区域平均
latmin=1;
latmax=25;
lonmin=101;
lonmax=120;

Nlat=find(LAT>=latmin&LAT<=latmax);
Nlon=find(LON>=lonmin&LON<=lonmax);

LON_SCS=LON(Nlon);
LAT_SCS=LAT(Nlat);
SSS_aer_SCS=SSS_aer(Nlon,Nlat,(1920-1920)+1:end);
SSS_ghg_SCS=SSS_hgh(Nlon,Nlat,(1920-1920)+1:end);
SSS_bmb_SCS=SSS_bmb(Nlon,Nlat,(1920-1920)+1:end);

SSS_SCS_ghg_a=SSS_ghg_SCS-nanmean(SSS_ghg_SCS,3);
SSS_SCS_aer_a=SSS_aer_SCS-nanmean(SSS_aer_SCS,3);
SSS_SCS_bmb_a=SSS_bmb_SCS-nanmean(SSS_bmb_SCS,3);
%区域平均
[nxx nyy ntt]=size(SSS_SCS_ghg_a);
tt=reshape(SSS_SCS_ghg_a,[nxx*nyy,ntt]);
SSS_SCS_ghg_a_time=squeeze(nanmean(tt));

[nxx nyy ntt]=size(SSS_SCS_aer_a);
tt=reshape(SSS_SCS_aer_a,[nxx*nyy,ntt]);
SSS_SCS_aer_a_time=squeeze(nanmean(tt));

[nxx nyy ntt]=size(SSS_SCS_bmb_a);
tt=reshape(SSS_SCS_bmb_a,[nxx*nyy,ntt]);
SSS_SCS_bmb_a_time=squeeze(nanmean(tt));
%trend
SSS_SCS_aer_trend=SSS_SCS_aer_a_time-detrend(SSS_SCS_aer_a_time);
SSS_SCS_ghg_trend=SSS_SCS_ghg_a_time-detrend(SSS_SCS_ghg_a_time);
SSS_SCS_bmb_trend=SSS_SCS_bmb_a_time-detrend(SSS_SCS_bmb_a_time);

[nyy ntt]=size(SSS_SCS_aer_a_time);
t2=[1:1:ntt];
y=squeeze(SSS_SCS_aer_a_time);
P=polyfit(t2',y',1);
SSS_SCS_aer_trend2=P(1);

[nyy ntt]=size(SSS_SCS_ghg_a_time);
t2=[1:1:ntt];
y=squeeze(SSS_SCS_ghg_a_time);
P=polyfit(t2',y',1);
SSS_SCS_ghg_trend2=P(1);

[nyy ntt]=size(SSS_SCS_bmb_a_time);
t2=[1:1:ntt];
y=squeeze(SSS_SCS_bmb_a_time);
P=polyfit(t2',y',1);
SSS_SCS_bmb_trend2=P(1);
%--------------------------------------------------------------------------
load('H:\SCS_SSS\图1SCS_SSS_trend_map\CESM_his_SCS_0_55m_trend.mat')

figure('Units','centimeters','Position',[8,2,21,8]);
plot(SSS_SCS_modelmean,'Color',[181/255,212/255,101/255],'LineWidth',2)
hold on
plot(SSS_SCS_ghg_a_time,'Color',[0.07,0.62,1.00],'LineWidth',2)
plot(SSS_SCS_aer_a_time,'Color',[1,99/255,71/255],'LineWidth',2)
plot(SSS_SCS_bmb_a_time,'Color',[147/255,112/255,219/255],'LineWidth',2)
plot([1 86],[0 0],'LineWidth',2,'Color',[.5 .5 .5])
plot(SSS_SCS_trend,'--','Color',[181/255,212/255,101/255],'LineWidth',2)
plot(SSS_SCS_aer_trend,'--','Color',[1,99/255,71/255],'LineWidth',2)
plot(SSS_SCS_ghg_trend,'--','Color',[0.07,0.62,1.00],'LineWidth',2)
plot(SSS_SCS_bmb_trend,'--','Color',[147/255,112/255,219/255],'LineWidth',2)
% plot(SSS_SCS_bmb_trend,'--g','LineWidth',1.5)
% grid on
set(gca,'Linewidth',1,'XTickLabel',...
    [1920 1930 1940 1950 1960 1970 1980 1990 2000],...
'XTick',1:10:86,'XLim',[1 86], 'yticklabel',-0.6:0.2:0.6,...
'YTick',-0.6:0.2:0.6, 'ylim',[-0.4 0.4],'Ycolor','k',...
'fontname','Times New Roman','fontsize',12,'fontweight','normal');
% text(2,-0.3,'CESM-GHGs Trend=-0.033 g/kg/100yr','fontsize',12)
% text(50,-0.3,'CESM-Aer Trend=0.41 g/kg/100yr','fontsize',12)
% text(50,-0.5,'CESM-bmb Trend=0.0017 g/kg/yr','fontsize',12)
hl=legend('CESM-LE(Trend=0.51)','GHGs(Trend=-0.033)','Aer(Trend=0.41)',...
    'bmb(Trend=0.17)');
set(hl,'box','off','fontsize',12)%'Orientation','horizon',
% print(gcf,'-dpng','-r600','CESM_his_SCS_0_55m_trend(time_line)_bmb')
% print(gcf,'-depsc2','-r600','CESM_S_0_55_single(time_line)')
% print -depsc2 -painters CESM_S_0_55_single(time_line)_bmb2.eps