clear all
%Various info 
dlon = 360;
dlat = 180;

%main pathes
O_PATH = '/cluster/shared/noresm/norcpm/Obs/SSH/CLS/';
FREE_PATH = '/cluster/home/ywang/NorESM1-CMIP6/noresm1-cmip6_historical_19500101/Ensemble_mean/ocn/';
ASSIM_PATH = '/cluster/home/ywang/NorCPM/norcpm-cmip6_analysis_19500115/Processed/ocn/';
ASSIM_PATH0 = '/cluster/home/ywang/NorCPM/noresm1-cmip6_analysis_19500115/Processed/ocn/';
Ave_PATH = '/cluster/shared/noresm/norcpm/Input/NorESM/f19_g16/CLS/';

%reading pivot
ipiv = ncgetvar([Ave_PATH 'pivots_SSH.nc'], 'ipiv');
jpiv = ncgetvar([Ave_PATH 'pivots_SSH.nc'], 'jpiv');
for j = 1:dlat
   lat(j) = -90 + j;
end
for i = 1:dlon
   lon(i) = (i - 1);
end

%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);
paream = mean(parea(:));
pareas = sum(sum(parea(:)));
plon=repmat(lon', 1, dlat);
plat=repmat(lat, dlon, 1);

if (~exist('Correlation_ssh_yearly.mat','file'))
cnt=0;
%reading data
for yr = 1993 : 2018
   yr
   cnt=cnt+1;
   %Data
   d_ssh(cnt,1:dlon,1:dlat) = ncread([O_PATH num2str(yr) '.nc'],'zo');

   %reading Free
   f_ssh = ncread([FREE_PATH 'free-average' num2str(yr) '.nc'],'sealv');
      
   %reading Assim
   a_ssh = ncread([ASSIM_PATH 'analysis-average' num2str(yr) '.nc'], 'sealv');

   %reading Assim
   a_ssh0 = ncgetvar([ASSIM_PATH0 'analysis-average' num2str(yr) '.nc'], 'sealv');

   %interp_model to obs
   f_ossh = zeros(dlon,dlat);
   nb_fossh = zeros(dlon,dlat);
   a_ossh = zeros(dlon,dlat);
   nb_aossh = zeros(dlon,dlat);
   a_ossh0 = zeros(dlon,dlat);
   nb_aossh0 = zeros(dlon,dlat);
   for i = 1:dlon
      for j = 1:dlat
         if ipiv(i,j) ~=0 & jpiv(i,j) ~=0
            f_ossh(i,j) = f_ossh(i,j) + f_ssh(ipiv(i,j),jpiv(i,j));
            nb_fossh(i,j) = nb_fossh(i,j) + 1;
            a_ossh(i,j) = a_ossh(i,j) + a_ssh(ipiv(i,j),jpiv(i,j));
            nb_aossh(i,j) = nb_aossh(i,j) + 1;
            a_ossh0(i,j) = a_ossh0(i,j) + a_ssh0(ipiv(i,j),jpiv(i,j));
            nb_aossh0(i,j) = nb_aossh0(i,j) + 1;
         end
      end
   end
   yr_f_ossh(cnt,1:dlon,1:dlat) = f_ossh ./ nb_fossh;
   yr_a_ossh(cnt,1:dlon,1:dlat) = a_ossh ./ nb_aossh;
   yr_a_ossh0(cnt,1:dlon,1:dlat) = a_ossh0 ./ nb_aossh0;
   yr_a_ossh3(cnt,1:dlon,1:dlat) = (yr_a_ossh(cnt,1:dlon,1:dlat) + yr_a_ossh0(cnt,1:dlon,1:dlat))/2;
end
%bias
bias_f=-squeeze(nanmean(d_ssh,1)-nanmean(yr_f_ossh,1));
bias_a=-squeeze(nanmean(d_ssh,1)-nanmean(yr_a_ossh,1));
bias_a0=-squeeze(nanmean(d_ssh,1)-nanmean(yr_a_ossh0,1));
bias_a3=-squeeze(nanmean(d_ssh,1)-nanmean(yr_a_ossh3,1));

for i = 1:dlon
   for j = 1:dlat
     %rmse_f(i,j) = 100*sqrt(nanmean(((yr_f_ossh(:,i,j)) - (d_ssh(:,i,j))).^2,1));
     %rmse_a(i,j) = 100*sqrt(nanmean(((yr_a_ossh(:,i,j)) - (d_ssh(:,i,j))).^2,1));
     %rmse_a0(i,j) = 100*sqrt(nanmean(((yr_a_ossh0(:,i,j)) - (d_ssh(:,i,j))).^2,1));
     corr_f(i,j) = corr((yr_f_ossh(:,i,j)), (d_ssh(:,i,j)));
     corr_a(i,j) = corr((yr_a_ossh(:,i,j)), (d_ssh(:,i,j)));
     corr_a0(i,j) = corr((yr_a_ossh0(:,i,j)), (d_ssh(:,i,j)));
     corr_a3(i,j) = corr((yr_a_ossh3(:,i,j)), (d_ssh(:,i,j)));
   end
end
   save(['Correlation_ssh_yearly.mat'],'corr_f','corr_a','corr_a0','corr_a3','bias_f','bias_a','bias_a0','bias_a3');
else
   load('Correlation_ssh_yearly.mat')
end
%rmse_f = 100*squeeze(sqrt(nanmean((yr_f_ossh - d_ssh).^2,1)/cnt));
%rmse_a = 100*squeeze(sqrt(nanmean((yr_a_ossh - d_ssh).^2,1)/cnt));
%rmse_a0= 100*squeeze(sqrt(nanmean((yr_a_ossh0 - d_ssh).^2,1)/cnt));
%bias_f = 100*squeeze(nanmean((yr_f_ossh - d_ssh),1));
%bias_a = 100*squeeze(nanmean((yr_a_ossh - d_ssh),1));
%bias_a0= 100*squeeze(nanmean((yr_a_ossh0 - d_ssh),1));

%%%%%%%%%%%%
%%%%Plot%%%%
%%%%%%%%%%%%
close all
if 1
figure(9)
clf
set(gcf, 'Renderer', 'opengl')
set(gcf, 'InvertHardCopy', 'off');
set(gca,'xticklabels',[],'yticklabels',[]);
%set(gca,'XAxisLocation','top','fontsize',16);
whitebg('w');
hold on
m_proj('Equidistant Cylindrical','lon',[-325 35],'lat',[-90 90]);
P=m_pcolor(plon, plat, bias_a);
set(P,'LineStyle','none')
P=m_pcolor(plon-360,plat, bias_a);
set(P,'LineStyle','none')
%title('dcppA-assim-i1','FontSize', 20);
caxis([-1 1]);
colormap('anom(20)');
m_grid;
m_ungrid;
m_coast('patch',[0.5 0.5 0.5]);
colorbar('h','FontSize', 16)
print('-depsc','bias_i1_SSH.eps')
end

if 1
figure(11)
clf
set(gcf, 'Renderer', 'opengl')
set(gcf, 'InvertHardCopy', 'off');
set(gca,'xticklabels',[],'yticklabels',[]);
%set(gca,'XAxisLocation','top','fontsize',16);
whitebg('w');
hold on
m_proj('Equidistant Cylindrical','lon',[-325 35],'lat',[-90 90]);
P=m_pcolor(plon, plat, abs(bias_a)-abs(bias_f));
set(P,'LineStyle','none')
P=m_pcolor(plon-360,plat, abs(bias_a)-abs(bias_f));
set(P,'LineStyle','none')
%title('dcppA-assim-i1','FontSize', 20);
caxis([-0.1 0.1]);
colormap('anom(20)');
m_grid;
m_ungrid;
m_coast('patch',[0.5 0.5 0.5]);
colorbar('h','FontSize', 16)
hold off
print('-depsc','bias_i1-free_SSH.eps')

figure(12)
clf
set(gcf, 'Renderer', 'opengl')
set(gcf, 'InvertHardCopy', 'off');
set(gca,'xticklabels',[],'yticklabels',[]);
%set(gca,'XAxisLocation','top','fontsize',16);
whitebg('w');
hold on
m_proj('Equidistant Cylindrical','lon',[-325 35],'lat',[-90 90]);
P=m_pcolor(plon, plat, abs(bias_a0)-abs(bias_a));
set(P,'LineStyle','none')
P=m_pcolor(plon-360,plat, abs(bias_a0)-abs(bias_a));
set(P,'LineStyle','none')
%title('dcppA-assim-i1','FontSize', 20);
caxis([-0.1 0.1]);
colormap('anom(20)');
m_grid;
m_ungrid;
m_coast('patch',[0.5 0.5 0.5]);
colorbar('h','FontSize', 16)
hold off
print('-depsc','bias_i2-i1_SSH.eps')
end