clc;clear all;close all
addpath('/data3/yiwen/cesm_final_exp/matlab_toolbox/m_map')
addpath('/data3/yiwen/cesm_final_exp/matlab_toolbox/colorbrewer')
pathi='/data3/yiwen/cesm_final_exp/plot_test/OBS_data/HadISST_sst.nc'

sst=ncread(pathi,'sst');
lon=ncread(pathi,'longitude');lat=ncread(pathi,'latitude');
time=ncread(pathi,'time');

% lon_loc=lon>=100 &lon<=156;lat_loc=lat>=10 &lat<=43;
lon_loc=lon>=100 &lon<=140;lat_loc=lat>=10 &lat<=40;

%1870-2021
for it = 1:1812/12
SST_yr(:,:,it)=nanmean(sst(:,:,(it-1)*12+1:(it-1)*12+12),3);
end

SST_yr_region=SST_yr(lon_loc,lat_loc,115:145);

time=linspace(1984,2014,31);
time2=linspace(1984,2000,17);
time3=linspace(2001,2014,14);

%1984-2014%
SST_avg1=nanmean(SST_yr_region,3);

SST_yr_2000=SST_yr_region(:,:,1:17); 
SST_avg2=nanmean(SST_yr_2000,3);

SST_yr_2015=SST_yr_region(:,:,18:31); 
SST_avg3=nanmean(SST_yr_2015,3);


%1994-1998
% SST_avg1=nanmean(SST_yr(lon_loc,lat_loc,124:128),3);

for i=1:40
  for j=1:30
     b1 = polyfit(time,squeeze(SST_yr_region(i,j,:))',1);
     A1(i,j)=b1(1);
     b2 = polyfit(time2,squeeze(SST_yr_2000(i,j,:))',1);
     A2(i,j)=b2(1);
     b3 = polyfit(time3,squeeze(SST_yr_2015(i,j,:))',1);
     A3(i,j)=b3(1);
  end
end

%linear trend :year mean
SST_linear=squeeze(nanmean(nanmean(SST_yr_region,1),2));
time=linspace(1984,2014,31);

b=polyfit(time,SST_linear',1);
fr=polyval(b,time);

plot(time,SST_linear','lineWidth',3);hold on
plot(time,fr,'lineWidth',3);
h1 = legend('Data', 'Linear trend', 'Location','NE')
[xlim ylim]
textposx = diff(xlim)*0.50+min(xlim);
textposy = diff(ylim)*0.95+min(ylim);
text(textposx,textposy, sprintf('y = %.3f*x %c %0.2f', b(1),...
     char(45-(sign(b(1))+1)), abs(b(2))), 'FontSize',20)
xlabel('year');ylabel('SST');title('Hadl SST trend')
set(gca,'Fontsize',20);set(gca,'Linewidth',3);

%%
[xx,yy]=meshgrid(lon(lon_loc),lat(lat_loc));

set(0,'defaultfigurecolor','w')
subplot(1,2,1)
m_proj('miller','lat',[10,40],'lon',[100,140]);
m_pcolor(xx,yy,SST_avg1');shading interp;hold on
m_coast('patch','w','edgecolor','k');
m_grid('linestyle','none','Fontsize',20,'Linewidth',3);
title({'SST average 1984-2014)'})
colormap(flipud(brewermap(30, 'RdYlBu')));
c = colorbar('location','southoutside');c.Title.String = "C";
set(c,'FontSize',15);caxis([15 30])
%colormap('jet')
set(gca,'Fontsize',20);set(gca,'Linewidth',3);

subplot(1,2,2)
m_proj('miller','lat',[10,40],'lon',[100,140]);
m_pcolor(xx,yy,A1');shading interp;hold on
m_coast('patch','w','edgecolor','k');
m_grid('linestyle','none','Fontsize',20,'Linewidth',3);
colormap(flipud(brewermap(30, 'RdYlBu')));
c = colorbar('location','southoutside');c.Title.String = "yearly trend";
set(c,'FontSize',15);caxis([-0.1 0.1])
title({'SST trend 1984-2014)'})
set(gca,'Fontsize',20);set(gca,'Linewidth',3);



