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

%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/';

for j = 1:dlat
   lat(j) = -90 + j;
end
for i = 1:dlon
   lon(i) = (i - 1);
end

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(1);clf;hold on
fld=bias_a;
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_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
m_proj('Equidistant Cylindrical','lon',[-325 35],'lat',[-90 90]);
P=m_pcolor(lon, lat, fld);
set(P,'LineStyle','none')
caxis([-1 1]);
colormap('anom(20)');
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
set(gca,'DataAspectRatio',[1.2 1 1e+16])
m_coast('patch',[0.5 0.5 0.5]);
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_i1_SSH.eps')
end




if 1
fld=abs(bias_a)-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
m_proj('Equidistant Cylindrical','lon',[-325 35],'lat',[-90 90]);
P=m_pcolor(lon, lat, fld);
set(P,'LineStyle','none')
caxis([-0.1 0.1]);
colormap('anom(20)');
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
set(gca,'DataAspectRatio',[1.2 1 1e+16])
m_coast('patch',[0.5 0.5 0.5]);
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_i1-free_SSH.eps')
end

if 1
fld=abs(bias_a0)-abs(bias_a);
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(3);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
m_proj('Equidistant Cylindrical','lon',[-325 35],'lat',[-90 90]);
P=m_pcolor(lon, lat, fld);
set(P,'LineStyle','none')
caxis([-0.1 0.1]);
colormap('anom(20)');
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
set(gca,'DataAspectRatio',[1.2 1 1e+16])
m_coast('patch',[0.5 0.5 0.5]);
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_i2-i1_SSH.eps')
end


if 0
figure(2)
clf
set(gcf, 'Renderer', 'opengl')
set(gcf, 'InvertHardCopy', 'off');
set(gca,'xticklabels',[],'yticklabels',[]);
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')
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('-depsc2','bias_i1-free_SSH.eps')

figure(2)
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('-depsc2','bias_i2-i1_SSH.eps')
end
