clear all
%Various info 
dlon=72;
dlat=36;
nx=320;
ny=384;
year_start=1950;
year_end=2018;
nb_yr=year_end-year_start + 1;
O_PATH='../../../../obsdata/regridded/SAT/';
FREE_PATH='/tos-project4/NS9039K/shared/norcpm/cases/NorESM1-CMIP6/noresm1-cmip6_historical_19500101/Ensemble_mean/atm/';
NorCPM_PATH1='/cluster/home/ywang/NorCPM/norcpm-cmip6_analysis_19500115/Processed/atm/';
NorCPM_PATH2='/cluster/home/ywang/NorCPM/noresm1-cmip6_analysis_19500115/Processed/atm/';

lat=ncread([O_PATH 'HadCRUT.4.6.0.0.median.nc'], 'latitude')';
lon=ncread([O_PATH 'HadCRUT.4.6.0.0.median.nc'], 'longitude')';
lon(lon<0)=lon(lon<0)+360;
%spatial distance
for j = 2:dlat-1
   dx = Haversin_dist(lon(1), lat(j), lon(3), lat(j)) / 2;
   dy = Haversin_dist(lon(1), lat(j-1), lon(1), lat(j+1)) / 2;
   parea(1:dlon, j) = dx * dy;
end
parea(1:dlon, 1) = parea(1:dlon, 2);
parea(1:dlon,dlat) = parea(1:dlon,dlat-1);
plon=repmat(lon', 1, dlat);
plat=repmat(lat, dlon, 1);

%calculate pivot
mlat= ncread([NorCPM_PATH1 'analysis-average2018-12.nc'], 'lat');
mlon= ncread([NorCPM_PATH1 'analysis-average2018-12.nc'], 'lon');

if exist(['T2m_yearly_NorCPM_HadCRUT4.mat'])
  load(['T2m_yearly_NorCPM_HadCRUT4.mat']);
else
  cnt = 1;
  %Data
  sst = squeeze(ncread([O_PATH 'HadCRUT.4.6.0.0.median.nc'],'temperature_anomaly'));
  clim=ncread([O_PATH 'absolute.nc'],'tem');
  %mask=find(sst < 0.);
  %sst(mask)=nan;
  for yr = year_start : year_end
   yr
   %Data
   ssto(cnt,:,:) = squeeze(nanmean(sst(:,:,1+(yr - 1850)*12:12+(yr - 1850)*12),3))+squeeze(nanmean(clim(:,:,:),3));
   %reading Free
   tmp = ncread([FREE_PATH 'free-average' num2str(yr) '.nc'],'TREFHT');
   sstf(cnt,:,:)=interp2(mlat',mlon,tmp,lat,lon');
   %reading Assim
   tmp = ncread([NorCPM_PATH1 'analysis-average' num2str(yr) '.nc'], 'TREFHT');
   sst1(cnt,:,:)=interp2(mlat',mlon,tmp,lat,lon');
   %reading Assim
   tmp = ncread([NorCPM_PATH2 'analysis-average' num2str(yr) '.nc'], 'TREFHT');
   sst2(cnt,:,:)=interp2(mlat',mlon,tmp,lat,lon');    
   %combination
   sst3(cnt,:,:)=(sst1(cnt,:,:)+sst2(cnt,:,:))/2;
   cnt = cnt + 1;
  end
  disp(cnt-1)
  bias_f=-squeeze(nanmean(ssto-sstf));
  bias_1=-squeeze(nanmean(ssto-sst1));  
  bias_2=-squeeze(nanmean(ssto-sst2));
  
  for i = 1:dlon
     for j = 1:dlat
         a=squeeze(sstf(:,i,j));
         b=squeeze(ssto(:,i,j));
         a=a(~isnan(b));
         b=b(~isnan(b));
         if ~isempty(b)
         R_f(i,j)=corr(a,b);
         else
         R_f(i,j)=nan;
         end
         a=squeeze(sst1(:,i,j));
         b=squeeze(ssto(:,i,j));
         a=a(~isnan(b));
         b=b(~isnan(b));
         if ~isempty(b)
         R_1(i,j)=corr(a,b);
         else
         R_1(i,j)=nan;
         end
         a=squeeze(sst2(:,i,j));
         b=squeeze(ssto(:,i,j));
         a=a(~isnan(b));
         b=b(~isnan(b));
         if ~isempty(b)
         R_2(i,j)=corr(a,b);
         else
         R_2(i,j)=nan;
         end
         a=squeeze(sst3(:,i,j));
         b=squeeze(ssto(:,i,j));
         a=a(~isnan(b));
         b=b(~isnan(b));
         if ~isempty(b)
         R_3(i,j)=corr(a,b);
         else
         R_3(i,j)=nan;
         end
         %rmse_f(i,j)=sqrt(nanmean((ssto(:,i,j)-sstf(:,i,j)).^2)-bias_f(i,j)^2);
         %rmse_1(i,j)=sqrt(nanmean((ssto(:,i,j)-sst1(:,i,j)).^2)-bias_1(i,j)^2);
         %rmse_2(i,j)=sqrt(nanmean((ssto(:,i,j)-sst2(:,i,j)).^2)-bias_2(i,j)^2);
     end
  end
  % save data
  save(['T2m_yearly_NorCPM_HadCRUT4.mat'],'R_f','R_1','R_2','bias_f','bias_1','bias_2');
end
%%%%%%%%%%%%
%%%%Plot%%%%
%%%%%%%%%%%%
close all
alon(1:36,:)=plon(37:end,:);
alon(37:72,:)=plon(1:36,:);
plon=alon;
c(1:36,:)=R_f(37:end,:);
c(37:72,:)=R_f(1:36,:);
R_f=c;
c(1:36,:)=R_1(37:end,:);
c(37:72,:)=R_1(1:36,:);
R_1=c;
c(1:36,:)=R_2(37:end,:);
c(37:72,:)=R_2(1:36,:);
R_2=c;
c(1:36,:)=bias_f(37:end,:);
c(37:72,:)=bias_f(1:36,:);
bias_f=c-273.15;
c(1:36,:)=bias_1(37:end,:);
c(37:72,:)=bias_1(1:36,:);
bias_1=c-273.15;
c(1:36,:)=bias_2(37:end,:);
c(37:72,:)=bias_2(1:36,:);
bias_2=c-273.15;

if 1
figure(4)
set(gcf, 'Renderer', 'opengl')
set(gcf, 'InvertHardCopy', 'off');
%set(gca,'xticklabels',[],'yticklabels',[]);
%set(gca,'XAxisLocation','top','fontsize',18);
whitebg('w');
hold on
m_proj('Equidistant Cylindrical','lon',[-325 35],'lat',[-90 90]);
m_coast('patch','w', 'edgecolor','k');
P=m_pcolor(plon, plat, bias_1);
set(P,'LineStyle','none')
P=m_pcolor(plon-360,plat, bias_1);
set(P,'LineStyle','none')
%title('Yearly SST RMSE','FontSize', 28);
caxis([-5 5]);
colormap(anom(20));
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
%m_coast('patch',[0.5 0.5 0.5]);
c=colorbar('h','FontSize', 20);
%xlabel(c,'K','FontSize', 22);
%print('-depsc','rmse_free_SSTA.eps')
print('-depsc','Bias_sat_i1.eps')
end

if 1
figure(6)
clf
set(gcf, 'Renderer', 'opengl')
set(gcf, 'InvertHardCopy', 'off');
%set(gca,'XAxisLocation','top','fontsize',18);
whitebg('w');
hold on
m_proj('Equidistant Cylindrical','lon',[-325 35],'lat',[-90 90]);
m_coast('patch','w', 'edgecolor','k');
P=m_pcolor(plon, plat, abs(bias_1)-abs(bias_f));
set(P,'LineStyle','none')
P=m_pcolor(plon-360,plat, abs(bias_1)-abs(bias_f));
set(P,'LineStyle','none')
%title('Yearly SST RMSE','FontSize', 28);
caxis([-1 1]);
colormap(anom(20));
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
%m_ungrid;
%m_coast('patch',[0.5 0.5 0.5]);
c=colorbar('h','FontSize', 20);
%xlabel(c,'K','FontSize', 22);
%print('-depsc','rmse_free_SSTA.eps')
print('-depsc','Bias_sat_i1-free.eps')

figure(7)
clf
set(gcf, 'Renderer', 'opengl')
set(gcf, 'InvertHardCopy', 'off');
%set(gca,'XAxisLocation','top','fontsize',18);
whitebg('w');
hold on
m_proj('Equidistant Cylindrical','lon',[-325 35],'lat',[-90 90]);
m_coast('patch','w', 'edgecolor','k');
P=m_pcolor(plon, plat, abs(bias_2)-abs(bias_1));
set(P,'LineStyle','none')
P=m_pcolor(plon-360,plat, abs(bias_2)-abs(bias_1));
set(P,'LineStyle','none')
%title('Yearly SST RMSE','FontSize', 28);
caxis([-1 1]);
colormap(anom(20));
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
%m_coast('patch',[0.5 0.5 0.5]);
c=colorbar('h','FontSize', 20);
%xlabel(c,'K','FontSize', 22);
%print('-depsc','rmse_free_SSTA.eps')
print('-depsc','Bias_sat_i2-i1.eps')
end