clear all
addpath ~/matlab/m_map
addpath ~/matlab/mats

%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;
plon=repmat(lon', 1, dlat);
plat=repmat(lat, dlon, 1);


if exist(['T2m_yearly_NorCPM_HadCRUT4.mat'])
  load(['T2m_yearly_NorCPM_HadCRUT4.mat']);
else
  mlat= ncread([NorCPM_PATH1 'analysis-average2018-12.nc'], 'lat');
  mlon= ncread([NorCPM_PATH1 'analysis-average2018-12.nc'], 'lon');
  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
fld=bias_1;
lon=plon-360;
lat=plat;
lon=cat(1,lon,lon+360);
lat=cat(1,lat,lat);
fld=cat(1,fld,fld);
lon=cat(2,lon(:,1),lon,lon(:,end));
lat=cat(2,lat(:,1)-2.5,lat,lat(:,end)+2.5);
fld=cat(2,fld(:,1),fld,fld(:,end));
figure(1);clf;hold on
figure_height_scale=0.65;
m_proj('Equidistant Cylindrical','lon',[-325 35],'lat',[-90 90]);
set(gcf,'paperunits','inches','papersize',[8 6*figure_height_scale], ...
        'paperposition',[0 0 8 6*figure_height_scale],'renderer','painters')
set(gca,'outerposition',[0 0 1 1],'position',[0.1 0.2 0.8 0.75])
set(gca,'outerposition',[0 0 1 1],'position',[0.1 0.15 0.8 0.85])
set(gcf, 'InvertHardCopy', 'off');
whitebg('w');
hold on
P=m_pcolor(lon, lat, fld);
set(P,'LineStyle','none')
m_coast('color','k');
caxis([-6 6]);
colormap(anom(20));
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
set(gca,'DataAspectRatio',[1.2 1 1e+16])
h=colorbar('southoutside','FontSize', 20);
%shrink=0.8;                                 
shrink=0.75;                                 
pos=get(h,'pos');                     
%set(h,'pos',[pos(1)+(1-shrink)/2*pos(3)  0.45*pos(2) shrink*pos(3) shrink*pos(4)]);
set(h,'pos',[pos(1)+(1-shrink)/2*pos(3)  0.4*pos(2) shrink*pos(3) shrink*pos(4)]);
print('-depsc2','Bias_sat_i1.eps')
end

if 1
fld=abs(bias_1)-abs(bias_f);
lon=plon-360;
lat=plat;
lon=cat(1,lon,lon+360);
lat=cat(1,lat,lat);
fld=cat(1,fld,fld);
lon=cat(2,lon(:,1),lon,lon(:,end));
lat=cat(2,lat(:,1)-2.5,lat,lat(:,end)+2.5);
fld=cat(2,fld(:,1),fld,fld(:,end));
figure(2);clf;hold on
figure_height_scale=0.65;
m_proj('Equidistant Cylindrical','lon',[-325 35],'lat',[-90 90]);
set(gcf,'paperunits','inches','papersize',[8 6*figure_height_scale], ...
        'paperposition',[0 0 8 6*figure_height_scale],'renderer','painters')
set(gca,'outerposition',[0 0 1 1],'position',[0.1 0.2 0.8 0.75])
set(gca,'outerposition',[0 0 1 1],'position',[0.1 0.15 0.8 0.85])
set(gcf, 'InvertHardCopy', 'off');
whitebg('w');
hold on
P=m_pcolor(lon, lat, fld);
set(P,'LineStyle','none')
m_coast('color','k');
caxis([-1 1]);
colormap(anom(20));
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
set(gca,'DataAspectRatio',[1.2 1 1e+16])
h=colorbar('southoutside','FontSize', 20);
%shrink=0.8;
shrink=0.75;
pos=get(h,'pos');
%set(h,'pos',[pos(1)+(1-shrink)/2*pos(3)  0.45*pos(2) shrink*pos(3) shrink*pos(4)]);
set(h,'pos',[pos(1)+(1-shrink)/2*pos(3)  0.4*pos(2) shrink*pos(3) shrink*pos(4)]);
print('-depsc2','Bias_sat_i1-free.eps')
end

if 1
fld=abs(bias_2)-abs(bias_1);
lon=plon-360;
lat=plat;
lon=cat(1,lon,lon+360);
lat=cat(1,lat,lat);
fld=cat(1,fld,fld);
lon=cat(2,lon(:,1),lon,lon(:,end));
lat=cat(2,lat(:,1)-2.5,lat,lat(:,end)+2.5);
fld=cat(2,fld(:,1),fld,fld(:,end));
figure(2);clf;hold on
figure_height_scale=0.65;
m_proj('Equidistant Cylindrical','lon',[-325 35],'lat',[-90 90]);
set(gcf,'paperunits','inches','papersize',[8 6*figure_height_scale], ...
        'paperposition',[0 0 8 6*figure_height_scale],'renderer','painters')
set(gca,'outerposition',[0 0 1 1],'position',[0.1 0.2 0.8 0.75])
set(gca,'outerposition',[0 0 1 1],'position',[0.1 0.15 0.8 0.85])
set(gcf, 'InvertHardCopy', 'off');
whitebg('w');
hold on
P=m_pcolor(lon, lat, fld);
set(P,'LineStyle','none')
m_coast('color','k');
caxis([-1 1]);
colormap(anom(20));
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
set(gca,'DataAspectRatio',[1.2 1 1e+16])
h=colorbar('southoutside','FontSize', 20);
%shrink=0.8;
shrink=0.75;
pos=get(h,'pos');
%set(h,'pos',[pos(1)+(1-shrink)/2*pos(3)  0.45*pos(2) shrink*pos(3) shrink*pos(4)]);
set(h,'pos',[pos(1)+(1-shrink)/2*pos(3)  0.4*pos(2) shrink*pos(3) shrink*pos(4)]);
print('-depsc2','Bias_sat_i2-i1.eps')
end


if 0
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

