clc
clear
close all;
ncdisp H:\data\CESM_his\PRECC\b.e11.B20TRC5CNBDRD.f09_g16.001.cam.h0.PRECC.185001-200512.nc
datadir='H:\data\CESM_his\PRECC\';
fielist=dir([datadir,'*.nc']);
k=length(fielist);
varname ='PRECC';
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_pre=lonData(Nlon);
    LAT_pre=latData(Nlat);
    
    stlo = min(Nlon); stla = min(Nlat);stti = 1;
    start = [stlo,stla,stti]; %创建start变量，每一维开始的变量
    locount = length(Nlon); lacount =length(Nlat); ticount = 1032;
    count = [locount,lacount,ticount]; %创建count变量，从每一维的start开始读取的总数目
    stride = [1,1,1]; %设置读取的步长
    PRECC= double(ncread(filename,varname,start,count,stride));
    PRECC_r(:,:,:,r)=PRECC*86400000;
    
%     PRECL_His_mean(:,:,r)=squeeze(nanmean(PRECL,3));
end
% save CESM_PRECC_SCS_his_trend.mat PRECC_r LON LAT
load('H:\SCS_SSS\图2SCS_PRE_trend_map\CESM_PRECC_SCS_his_trend.mat')
 load('H:\SCS_SSS\图2SCS_PRE_trend_map\CESM_PRECL_SCS_his_trend.mat')
%trend
PRE_r=PRECC_r+PRECL_r;
[nxx nyy ntt nrr]=size(PRE_r);
%annual 
for i=1:ntt./12
   PRE_an(:,:,i,:)=squeeze(nanmean(PRE_r(:,:,(i-1)*12+1:(i-1)*12+12,:),3));
end
PRE_an_1=PRE_an(:,:,(1920-1920)+1:(2005-1920)+1,:);
[nyy nzz ntt nrr]=size(PRE_an_1);
t2=[1:1:ntt];
for n=1:nrr
    for i=1:nyy
        for j=1:nzz
            y=squeeze(PRE_an_1(i,j,:,n));
            P=polyfit(t2,y',1);
            PRE_his_trend_r(i,j,n)=P(1);
        end
    end
end
%modelmean
PRE_his_trend=squeeze(nanmean(PRE_his_trend_r,3));

%信度检验
[nxx nyy nrr]=size(PRE_his_trend_r);
for i=1:nxx
    for j=1:nyy
        PRE_trend_tt(i,j)=ttest(PRE_his_trend_r(i,j,:),0,'Alpha',0.01);
    end
end
PRE_trend=squeeze(nanmean(PRE_his_trend_r,3));

PRE_trend_new=PRE_trend;
PRE_trend_new(PRE_trend_tt==0)=NaN;

figure('Units','centimeters','Position',[8,2,21,15]);
m_proj('Mercator','lat',[min(LAT_pre) max(LAT_pre)],'lon',[min(LON_pre) max(LON_pre)]);
[long,latt]=meshgrid(double(LON_pre),double(LAT_pre));
m_contourf(long,latt,PRE_his_trend'*100,500,'linestyle','none');
m_coast('color','[.5 .5 .5]','linewidth',1.5);
% colorbar;
hold on
[nxx nyy]=size(PRE_trend_new);
for i=1:nxx
    for j=1:nyy
         if(~isnan(PRE_trend_new(i,j)))
             m_plot(long(j,i),latt(j,i),'k.','markersize',3.5);
         end
    end
end
% 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)
m_grid('xtick',10,'ytick',7,'tickdir','out','yaxislocation','left','fontsize',10);%grid设置
caxis([-2.9,2.9]);
title('(a) CESM-LE Precip trend 1920-2005','fontsize',15)


load('J:/盐度/年代际变化/color/gran_blue_ncl.mat')
L=length(colors);
L2=floor(L/2);
a=colors(1:L2-10,:);b=colors(L2+10:end,:);
c=[a;b];
colormap(c);
hcb=colorbar('location','SouthOutside','position',[0.1,0.12,0.85,0.025]);
%  print(gcf,'-dpng','-r600','CESM_PREC_SCS_trend(1920_1980)');
%  print(gcf,'-depsc2','-r600','CESM_PREC_SCS_trend(1920_2005)')
% print -depsc2 -painters CESM_PREC_SCS_trend(1920_2005).eps
%%
%--------------------------------------------------------------------------
latmin=1;
latmax=25;
lonmin=101;
lonmax=120;

Nlat=find(LAT_pre>=latmin&LAT_pre<=latmax);
Nlon=find(LON_pre>=lonmin&LON_pre<=lonmax);

LON_SCS_PRE=LON_pre(Nlon);
LAT_SCS_PRE=LAT_pre(Nlat);

prec_an_SCS=PRE_an_1(Nlon,Nlat,:,:);

prec_an_SCS_a=prec_an_SCS-nanmean(prec_an_SCS,3);

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

%modelmean
prec_SCS_modelmean=squeeze(nanmean(prec_an_SCS_a_time,2));
%trend
prec_SCS_trend=prec_SCS_modelmean-detrend(prec_SCS_modelmean);

[nyy ntt]=size(prec_SCS_modelmean);
t2=[1:1:nyy];
y=squeeze(prec_SCS_modelmean);
P=polyfit(t2',y',1);
prec_an_SCS_trend_2=P(1);
%-------mask-------------------------------------------------------------------
load('I:\IO_Dome\CESM_his\SALT\His_SCS\SCS_mark.mat')

[XX,YY]=meshgrid(LAT_SCS,LON_SCS);
[X,Y]=meshgrid(LAT_SCS_PRE,LON_SCS_PRE);%插值之后的LON_E,LAT_E
precip_mask=griddata(XX,YY,SCS_mark,X,Y);

PP_SCS_new=prec_an_SCS;
[nxx nyy ntt nrr]=size(prec_an_SCS);
precip_mask_new=repmat(precip_mask,[1,1,ntt,nrr]);
PP_SCS_new(isnan(precip_mask_new))=NaN;

PP_SCS_new_a=PP_SCS_new-nanmean(PP_SCS_new,3);  
%区域平均
[nxx nyy ntt nrr]=size(PP_SCS_new_a);
tt=reshape(PP_SCS_new_a,[nxx*nyy,ntt,nrr]);
PP_SCS_new_a_time=squeeze(nanmean(tt));

%modelmean
PP_new_SCS_modelmean=squeeze(nanmean(PP_SCS_new_a_time,2));

PP_SCS_new_trend=PP_new_SCS_modelmean-detrend(PP_new_SCS_modelmean);

[nyy ntt]=size(PP_new_SCS_modelmean);
t2=[1:1:nyy];
y=squeeze(PP_new_SCS_modelmean);
P=polyfit(t2',y',1);
PP_SCS_trend_2=P(1);
%--------------------------------------------------------------------------
figure('Units','centimeters','Position',[8,2,21,8]);
plot(prec_SCS_modelmean,'Color',[1 99/255,71/255],'LineWidth',2)
hold on
plot(PP_new_SCS_modelmean,'Color',[0.07,0.62,1.00],'LineWidth',2)
plot(prec_SCS_trend,'--','LineWidth',2,'Color',[1 99/255,71/255])
plot(PP_SCS_new_trend,'--','LineWidth',2,'Color',[0.07,0.62,1.00])
% plot(PP_SCS_new_a_time,'Color', [0.5 0.5 0.5],'LineWidth',1.5)
plot([1 86],[0 0],'LineWidth',2,'Color',[.5 .5 .5])
plot(prec_SCS_modelmean,'Color',[1 99/255,71/255],'LineWidth',2)
plot(PP_new_SCS_modelmean,'Color',[0.07,0.62,1.00],'LineWidth',2)
plot(prec_SCS_trend,'--','LineWidth',2,'Color',[1 99/255,71/255])
plot(PP_SCS_new_trend,'--','LineWidth',2,'Color',[0.07,0.62,1.00])
% 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.8:0.2:0.8,...
'YTick',-0.8:0.2:0.8, 'ylim',[-0.8 0.8],'Ycolor','k',...
'fontname','Times New Roman','fontsize',12,'fontweight','normal');
% text(1.5,-0.5,'SCS Trend=-0.99 mm/day/100yr','fontsize',12)
% text(1.5,-0.7,'SCS Ocean Trend=-0.83 mm/day/100yr ','fontsize',12)
hl=legend('SCS(Trend=-0.99 mm/day/100yr)','SCS Ocean(Trend=-0.83 mm/day/100yr)');
set(hl,'box','off','fontsize',15) %'Orientation','horizon', 
text(2,0.7,'b','fontsize',20)
%  print(gcf,'-dpng','-r600','CESM_PREC_SCS_trend(time_line1920_2005)');
%  print(gcf,'-depsc2','-r600','CESM_PREC_SCS_trend(time_line1920_2005)')
% print -depsc2 -painters CESM_PREC_SCS_trend(time_line1920_2005).eps
%前后10年差
%前10年
PRE_q10=squeeze(nanmean(PP_SCS_new_a_time(1:10)))./55*35;

%前后差
PRE_c=squeeze(nanmean(PP_SCS_new_a_time(end-9:end)))./55*35-PRE_q10;