clc
clear
close all;
ncdisp E:\data\CESM_his\SLAT_global\SALT_global_annual\SALT_historical_global_annual_001.pop.nc
datadir='E:\data\CESM_his\SLAT_global\SALT_global_annual\';
fielist=dir([datadir,'*.nc']);
k=length(fielist);
varname ='SALT';

for r=1:k
    r
    filename=[datadir,fielist(r).name];
    
    %     sstData=double(ncread(filename,'TEMP'));
    lonData=double(ncread(filename,'lon'));
    latData=double(ncread(filename,'lat'));
    levData=double(ncread(filename,'level'))./100;
%     time=double(ncread(filename,'time'));
%     time=time+datenum(0000,1,1,0,0,0);%定位时间,原始时间为对应1800,1,1的天数
%     date=datestr(time,'yyyymmdd');%将时间转化为年月日的形式
    
    latmin=-10;
    latmax=30;
    lonmin=90;
    lonmax=180;
    
    Nlat=find(latData>=latmin&latData<=latmax);
    Nlon=find(lonData>=lonmin&lonData<=lonmax);
    LON=lonData(Nlon);
    LAT=latData(Nlat);
    
    stlo = min(Nlon); stla = min(Nlat);szzi=1;stti = 1;
    start = [stlo,stla,szzi,stti]; %创建start变量，每一维开始的变量
    locount = length(Nlon); lacount =length(Nlat);zicount=6; ticount = 86;
    count = [locount,lacount,zicount,ticount]; %创建count变量，从每一维的start开始读取的总数目
    stride = [1,1,1,1]; %设置读取的步长
    S= ncread(filename,varname,start,count,stride);
    S_0_55_r(:,:,:,r)=squeeze(nanmean(S,3)); 
    S_0_25_r(:,:,:,r)=squeeze(nanmean(S(:,:,1:3,:),3)); 
    S_55_r(:,:,:,r)=squeeze(nanmean(S(:,:,6,:),3));
end

%%
%表层
SSS=squeeze(S_0_55_r(:,:,(1920-1920)+1:(1980-1920)+1,:));
SSSa=SSS-nanmean(SSS,3);
%trend
[nyy nzz ntt nrr]=size(SSSa);
t2=[1:1:ntt];
for n=1:nrr
    for i=1:nyy
        for j=1:nzz
            y=squeeze(SSSa(i,j,:,n));
            P=polyfit(t2,y',1);
            S_trend_m(i,j,n)=P(1);
        end
    end
end
%信度检验
[nxx nyy nrr]=size(S_trend_m);
for i=1:nxx
    for j=1:nyy
        SSS_trend_tt(i,j)=ttest(S_trend_m(i,j,:),0,'Alpha',0.05);
    end
end
SSS_trend=squeeze(nanmean(S_trend_m,3));

SSS_trend_new=SSS_trend;
SSS_trend_new(SSS_trend_tt==0)=NaN;


% 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_SCS_ave=squeeze(nanmean(nanmean(SSS(Nlon,Nlat,:,:),3),4));
SSS_SCS_ave=squeeze(nanmean(nanmean(SSS,3),4));

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_trend'*100,500,'linestyle','none');
hold on;
m_coast('patch',[.7 .7 .7],'linewidth',1,'edgecolor',[.5 .5 .5])
[nxx nyy]=size(SSS_trend_new);
for i=1:nxx
    for j=1:nyy
         if(~isnan(SSS_trend_new(i,j)))
             m_plot(long(j,i),latt(j,i),'k.','markersize',3.5);
         end
    end
end
% [long,latt]=meshgrid(double(LON_SCS),double(LAT_SCS));
[c h]=m_contour(long,latt,SSS_SCS_ave',[28:0.5:35],'b','linewidth',1.5);
clabel(c,h,'LabelSpacing',700,'Color','b')
m_grid('xtick',10,'ytick',7,'tickdir','out','yaxislocation','left','fontsize',15);%grid设置
caxis([-1 1]);
lon_x1=101:120;lat_y1=ones(size(lon_x1));
m_plot(lon_x1,lat_y1,'k','linewidth',2)
lat_y2=1:0.1:25.5;lon_x2=ones(size(lat_y2))*101;
m_plot(lon_x2,lat_y2,'k','linewidth',2)
lat_y3=1: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)
%Green 120°E -127°E 15°N -25.5°N
lon_x5=120:127;lat_y5=ones(size(lon_x5))*25.5;
m_plot(lon_x5,lat_y5,'color',[0 139/255 0],'linewidth',3)
lon_x5=120:127;lat_y5=ones(size(lon_x5))*15;
m_plot(lon_x5,lat_y5,'color',[0 139/255 0],'linewidth',3)
lat_y6=15:26;lon_x6=ones(size(lat_y6))*127;
m_plot(lon_x6,lat_y6,'color',[0 139/255 0],'linewidth',3)
%Blue 103°E -109°E 3.5°S -1°N
lat_y6=-3.5:1;lon_x6=ones(size(lat_y6))*103;
m_plot(lon_x6,lat_y6,'color',[30/255 144/255 1],'linewidth',3)
lat_y6=-3.5:1;lon_x6=ones(size(lat_y6))*109;
m_plot(lon_x6,lat_y6,'color',[30/255 144/255 1],'linewidth',3)
lon_x5=103:109;lat_y5=ones(size(lon_x5))*-3.5;
m_plot(lon_x5,lat_y5,'color',[30/255 144/255 1],'linewidth',3)
% rectangle('Position',[-0.14,-0.11,0.4,0.54],'LineWidth',2,'LineStyle','-');
title('(a) CESM-his Salinity 0-55m 1920-1980 trend','fontsize',18)

load('E:\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('location','SouthOutside','position',[0.3,0.15,0.45,0.015]);
% print(gcf,'-dpng','-r600','CESM_his_SCS_SALT_trend(1920_1980)_1');
% print(gcf,'-depsc2','-r600','CESM_his_SCS_SALT_trend(1920_2005)')
print -depsc2 -painters CESM_his_SCS_SALT_trend(1920_1980).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_SCS=S_0_55_r(Nlon,Nlat,(1920-1920)+1:end,:);

SCS_mark=squeeze(nanmean(nanmean(SSS_SCS,3),4));
SCS_mark(~isnan(SCS_mark))=1;
% save SCS_mark.mat SCS_mark LON_SCS LAT_SCS

SSS_SCS_new=S_0_55_r(Nlon,Nlat,(1920-1920)+1:end,:);
[nxx nyy ntt nrr]=size(S_0_55_r(:,:,(1920-1920)+1:end,:));
SCS_mark_new=repmat(SCS_mark,[1,1,ntt,nrr]);
SSS_SCS_new(isnan(SCS_mark_new))=NaN;

SSS_SCS_a=SSS_SCS-nanmean(SSS_SCS,3);
SSS_SCS_new_a=SSS_SCS_new-nanmean(SSS_SCS_new,3);  

%区域平均
[nxx nyy ntt nrr]=size(SSS_SCS_a);
tt=reshape(SSS_SCS_a,[nxx*nyy,ntt,nrr]);
SSS_SCS_a_time=squeeze(nanmean(tt));

[nxx nyy ntt nrr]=size(SSS_SCS_new_a);
tt=reshape(SSS_SCS_new_a,[nxx*nyy,ntt,nrr]);
SSS_SCS_new_a_time=squeeze(nanmean(tt));
%modelmean
SSS_SCS_modelmean=squeeze(nanmean(SSS_SCS_a_time,2));
SSS_SCS_new_modelmean=squeeze(nanmean(SSS_SCS_new_a_time,2));
%trend
SSS_SCS_trend=SSS_SCS_modelmean-detrend(SSS_SCS_modelmean);
SSS_SCS_new_trend=SSS_SCS_new_modelmean-detrend(SSS_SCS_new_modelmean);

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

[nyy ntt]=size(SSS_SCS_new_modelmean);
t2=[1:1:nyy];
y=squeeze(SSS_SCS_new_modelmean');
P=polyfit(t2',y',1);
SSS_SCS_new_trend_2=P(1);
%--------------------------------------------------------------------------
figure('Units','centimeters','Position',[8,2,21,8]);
plot(SSS_SCS_modelmean,'k','LineWidth',1)
hold on
plot(SSS_SCS_a_time,'Color',[1 99/255,71/255],'LineWidth',2)%red
plot(SSS_SCS_new_modelmean,'b','LineWidth',1)
plot(SSS_SCS_trend,'--','Color',[1 99/255,71/255],'LineWidth',2)
plot(SSS_SCS_new_trend_2,'--b','LineWidth',1.5)
grid on
set(gca,'Linewidth',1,'XTickLabel',...
    [1940 1950 1960 1970 1980 1990 2000],...
'XTick',1:10:66,'XLim',[1 66], 'yticklabel',-0.6:0.2:0.6,...
'YTick',-0.6:0.2:0.6, 'ylim',[-0.6 0.6],'Ycolor','k',...
'fontname','Times New Roman','fontsize',12,'fontweight','normal');
text(20,-0.2,'SCS black Trend=0.0051 g/kg/yr','fontsize',12)
text(20,-0.4,'SCS mask Trend=0.0051 g/kg/yr','fontsize',12)

hl=legend('SCS black','SCS mask');
set(hl,'Orientation','horizon','box','off','fontsize',15)  
% print(gcf,'-dpng','-r600','CESM_his_SCS_0_55m_trend')

% save CESM_his_SCS_0_55m_trend.mat S_0_55_r S_55_r...
%     SSS_SCS_modelmean SSS_SCS_a_time SSS_SCS_trend SSS_SCS LON_SCS LAT_SCS
