%function Calculate_yearly_heat_salt_content_Anom_fanf(DEPTH_LIM) % Depth_lim equal to 200m or 500m for example
DEPTH_LIM=300;
%Various info 
dlon=360;
dlat=173;
nx=320
ny=384
year_start=1950
year_end=2018
O_PATH='/cluster/shared/noresm/norcpm/Obs/SST/HADISST2/'
TRUE_PATH='/cluster/shared/noresm/norcpm/Obs/ANALYSIS/EN.4.2.1/'
FREE_PATH='/tos-project4/NS9039K/shared/norcpm/cases/NorESM1-CMIP6/noresm1-cmip6_historical_19500101/Ensemble_mean/ocn/'
NorCPM_PATH1='/cluster/home/ywang/NorCPM/norcpm-cmip6_analysis_19500115/Processed/ocn/'
NorCPM_PATH2='/cluster/home/ywang/NorCPM/noresm1-cmip6_analysis_19500115/Processed/ocn/'
Input_PATH = '/cluster/shared/noresm/norcpm/Input/NorESM/f19_g16/EN421/'
grid_PATH='/cluster/shared/noresm/inputdata/ocn/micom/gx1v6/20101119/grid.nc'
ipiv=ncgetvar([Input_PATH 'pivots_EN4.nc'],'ipiv');
jpiv=ncgetvar([Input_PATH 'pivots_EN4.nc'],'jpiv');
for j=1:dlat
   lat(j)=-84+j;
end
for i=1:dlon
   lon(i)=i;
end
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,173)=parea(1:dlon,172);
paream=nanmean(parea(:));
plon=repmat(lon',1,dlat);
plat=repmat(lat,dlon,1);

if (~exist(['Heat_Salt_contents_yearly_' num2str(DEPTH_LIM) 'm_EN4_NorCPM.mat'],'file'))
   cnt=1
   hc_NorCPM1=zeros(dlon,dlat,(year_end-year_start+1));
   hc_NorCPM2=zeros(dlon,dlat,(year_end-year_start+1));
   hc_free=zeros(dlon,dlat,(year_end-year_start+1));
   hc_EN4=zeros(dlon,dlat,(year_end-year_start+1));
   sc_NorCPM1=zeros(dlon,dlat,(year_end-year_start+1));
   sc_NorCPM2=zeros(dlon,dlat,(year_end-year_start+1));
   sc_free=zeros(dlon,dlat,(year_end-year_start+1));
   sc_EN4=zeros(dlon,dlat,(year_end-year_start+1));

% compute globalmean.mean steric height 
weightsflg=0;
for yr=year_start:year_end
 yr 
    if weightsflg == 0 % Done just once because depths are always the same into these files
        weightsflg=1;
        depthbnds_NorCPM=ncgetvar([FREE_PATH 'free-average1980.nc'],'depth_bnds');
        kdm=size(depthbnds_NorCPM,2);
        % Reading the EN4 values
        depthbnds_EN4=ncgetvar([TRUE_PATH 'EN.4.2.1.f.analysis.g10.' num2str(yr) '.nc'],'depth_bnds');
        kkdm=size(depthbnds_EN4,2);
    end 
    
    temp_EN4=ncgetvar([TRUE_PATH 'EN.4.2.1.f.analysis.g10.' num2str(yr) '.nc'],'temperature')-273.15;    
%    avg_temp_EN4=ncgetvar([TRUE_PATH 'Anomaly/EN.4.1.1_1950-2009.nc'],'temperature');
%    temp_EN4=temp_EN4-avg_temp_EN4;
    saln_EN4=ncgetvar([TRUE_PATH 'EN.4.2.1.f.analysis.g10.' num2str(yr) '.nc'],'salinity');
%    avg_saln_EN4=ncgetvar([TRUE_PATH 'Anomaly/EN.4.1.1_1950-2009.nc'],'salinity');
%    saln_EN4=saln_EN4-avg_saln_EN4;

    fldt_NorCPM1=ncgetvar([NorCPM_PATH1 'analysis-average' num2str(yr) '.nc'],'templvl');
    fldt_NorCPM2=ncgetvar([NorCPM_PATH2 'analysis-average' num2str(yr) '.nc'],'templvl');
    fldt_free=ncgetvar([FREE_PATH 'free-average' num2str(yr) '.nc'],'templvl');
    flds_NorCPM1=ncgetvar([NorCPM_PATH1 'analysis-average' num2str(yr) '.nc'],'salnlvl');
    flds_NorCPM2=ncgetvar([NorCPM_PATH2 'analysis-average' num2str(yr) '.nc'],'salnlvl');
    flds_free=ncgetvar([FREE_PATH 'free-average' num2str(yr) '.nc'],'salnlvl');
    %Interpolate model to EN4
    temp_NorCPM1=zeros(dlon,dlat,kdm);
    temp_NorCPM2=zeros(dlon,dlat,kdm);
    temp_free=zeros(dlon,dlat,kdm);
    saln_NorCPM1=zeros(dlon,dlat,kdm);
    saln_NorCPM2=zeros(dlon,dlat,kdm);
    saln_free=zeros(dlon,dlat,kdm);
    nb_o=zeros(dlon,dlat,kdm);
    for i = 1:dlon
      for j = 1:dlat
        if ipiv(i,j) ~=0 & jpiv(i,j) ~=0
          temp_NorCPM1(i,j,:) = temp_NorCPM1(i,j,:) + fldt_NorCPM1(ipiv(i,j),jpiv(i,j),:);
          temp_NorCPM2(i,j,:) = temp_NorCPM2(i,j,:) + fldt_NorCPM2(ipiv(i,j),jpiv(i,j),:);
          temp_free(i,j,:) = temp_free(i,j,:) + fldt_free(ipiv(i,j),jpiv(i,j),:);
          saln_NorCPM1(i,j,:) = saln_NorCPM1(i,j,:) + flds_NorCPM1(ipiv(i,j),jpiv(i,j),:);
          saln_NorCPM2(i,j,:) = saln_NorCPM2(i,j,:) + flds_NorCPM2(ipiv(i,j),jpiv(i,j),:);
          saln_free(i,j,:) = saln_free(i,j,:) + flds_free(ipiv(i,j),jpiv(i,j),:);
          nb_o(i,j,:) = nb_o(i,j,:) + 1;
        end
      end
    end
    temp_NorCPM1 = temp_NorCPM1 ./ nb_o;
    temp_NorCPM2 = temp_NorCPM2 ./ nb_o;
    temp_free = temp_free ./ nb_o;
    saln_NorCPM1 = saln_NorCPM1 ./ nb_o;
    saln_NorCPM2 = saln_NorCPM2 ./ nb_o;
    saln_free = saln_free ./ nb_o;
    
    % Calculating of the heat/salt content
    for k=1:kdm
      if (depthbnds_NorCPM(2,k) <= DEPTH_LIM)
         hc_NorCPM1(:,:,cnt)=hc_NorCPM1(:,:,cnt)+temp_NorCPM1(:,:,k)*(depthbnds_NorCPM(2,k)-depthbnds_NorCPM(1,k));
         hc_NorCPM2(:,:,cnt)=hc_NorCPM2(:,:,cnt)+temp_NorCPM2(:,:,k)*(depthbnds_NorCPM(2,k)-depthbnds_NorCPM(1,k));
         hc_free(:,:,cnt)=hc_free(:,:,cnt)+temp_free(:,:,k)*(depthbnds_NorCPM(2,k)-depthbnds_NorCPM(1,k));
         sc_NorCPM1(:,:,cnt)=sc_NorCPM1(:,:,cnt)+saln_NorCPM1(:,:,k)*(depthbnds_NorCPM(2,k)-depthbnds_NorCPM(1,k));
         sc_NorCPM2(:,:,cnt)=sc_NorCPM2(:,:,cnt)+saln_NorCPM2(:,:,k)*(depthbnds_NorCPM(2,k)-depthbnds_NorCPM(1,k));
         sc_free(:,:,cnt)=sc_free(:,:,cnt)+saln_free(:,:,k)*(depthbnds_NorCPM(2,k)-depthbnds_NorCPM(1,k));
      else if (depthbnds_NorCPM(1,k) < DEPTH_LIM)
         hc_NorCPM1(:,:,cnt)=hc_NorCPM1(:,:,cnt)+temp_NorCPM1(:,:,k)*(DEPTH_LIM-depthbnds_NorCPM(1,k)); % We stop at DEPTH_LIM
         hc_NorCPM2(:,:,cnt)=hc_NorCPM2(:,:,cnt)+temp_NorCPM2(:,:,k)*(DEPTH_LIM-depthbnds_NorCPM(1,k)); % We stop at DEPTH_LIM
         hc_free(:,:,cnt)=hc_free(:,:,cnt)+temp_free(:,:,k)*(DEPTH_LIM-depthbnds_NorCPM(1,k)); % We stop at DEPTH_LIM
         sc_NorCPM1(:,:,cnt)=sc_NorCPM1(:,:,cnt)+saln_NorCPM1(:,:,k)*(DEPTH_LIM-depthbnds_NorCPM(1,k)); % We stop at DEPTH_LIM
         sc_NorCPM2(:,:,cnt)=sc_NorCPM2(:,:,cnt)+saln_NorCPM2(:,:,k)*(DEPTH_LIM-depthbnds_NorCPM(1,k)); % We stop at DEPTH_LIM
         sc_free(:,:,cnt)=sc_free(:,:,cnt)+saln_free(:,:,k)*(DEPTH_LIM-depthbnds_NorCPM(1,k)); % We stop at DEPTH_LIM
           end
      end
    end
    for kk3=1:kkdm
      if (depthbnds_EN4(2,kk3) <= DEPTH_LIM)
         hc_EN4(:,:,cnt)=hc_EN4(:,:,cnt)+temp_EN4(:,:,kk3)*(depthbnds_EN4(2,kk3)-depthbnds_EN4(1,kk3));
         sc_EN4(:,:,cnt)=sc_EN4(:,:,cnt)+saln_EN4(:,:,kk3)*(depthbnds_EN4(2,kk3)-depthbnds_EN4(1,kk3));
      else if (depthbnds_EN4(1,kk3) < DEPTH_LIM)
         hc_EN4(:,:,cnt)=hc_EN4(:,:,cnt)+temp_EN4(:,:,kk3)*(DEPTH_LIM-depthbnds_EN4(1,kk3));
         sc_EN4(:,:,cnt)=sc_EN4(:,:,cnt)+saln_EN4(:,:,kk3)*(DEPTH_LIM-depthbnds_EN4(1,kk3));
           end
      end
    end
    cnt=cnt+1;
end

save(['Heat_Salt_contents_yearly_' num2str(DEPTH_LIM) 'm_EN4_NorCPM.mat'],'hc_NorCPM1','hc_NorCPM2','hc_free','hc_EN4','sc_NorCPM1','sc_NorCPM2','sc_free','sc_EN4');

else
   load(['Heat_Salt_contents_yearly_' num2str(DEPTH_LIM) 'm_EN4_NorCPM.mat'])
end

if (~exist(['Correlations_hc_sc_yearly_' num2str(DEPTH_LIM) 'm_EN4_NorCPM.mat']))
   R_hc1=zeros(dlon,dlat);
   R_hc2=zeros(dlon,dlat);
   R_hcf=zeros(dlon,dlat);
   R_hc3=zeros(dlon,dlat);
   R_sc1=zeros(dlon,dlat);
   R_sc2=zeros(dlon,dlat);
   R_scf=zeros(dlon,dlat);
   R_sc3=zeros(dlon,dlat);
   rmse_hc1=zeros(dlon,dlat);
   rmse_hc2=zeros(dlon,dlat);
   rmse_hcf=zeros(dlon,dlat);
   rmse_sc1=zeros(dlon,dlat);
   rmse_sc2=zeros(dlon,dlat);
   rmse_scf=zeros(dlon,dlat);   
   bias_hc1=zeros(dlon,dlat);
   bias_hc2=zeros(dlon,dlat);
   bias_hcf=zeros(dlon,dlat);
   bias_sc1=zeros(dlon,dlat);
   bias_sc2=zeros(dlon,dlat);
   bias_scf=zeros(dlon,dlat);

   for m=1:dlon
      for n=1:dlat
       bias_hc1(m,n)=-nanmean(hc_EN4(m,n,:)-hc_NorCPM1(m,n,:))/DEPTH_LIM;
       bias_hc2(m,n)=-nanmean(hc_EN4(m,n,:)-hc_NorCPM2(m,n,:))/DEPTH_LIM;
       bias_hcf(m,n)=-nanmean(hc_EN4(m,n,:)-hc_free(m,n,:))/DEPTH_LIM;
       bias_sc1(m,n)=-nanmean(sc_EN4(m,n,:)-sc_NorCPM1(m,n,:))/DEPTH_LIM;
       bias_sc2(m,n)=-nanmean(sc_EN4(m,n,:)-sc_NorCPM2(m,n,:))/DEPTH_LIM;
       bias_scf(m,n)=-nanmean(sc_EN4(m,n,:)-sc_free(m,n,:))/DEPTH_LIM;  

       R_hc1(m,n)=corr(squeeze(hc_EN4(m,n,:)),squeeze(hc_NorCPM1(m,n,:)));
       R_hc2(m,n)=corr(squeeze(hc_EN4(m,n,:)),squeeze(hc_NorCPM2(m,n,:)));
       R_hcf(m,n)=corr(squeeze(hc_EN4(m,n,:)),squeeze(hc_free(m,n,:)));
       R_hc3(m,n)=corr(squeeze(hc_EN4(m,n,:)),squeeze(hc_NorCPM1(m,n,:)+hc_NorCPM2(m,n,:))/2);

       R_sc1(m,n)=corr(squeeze(sc_EN4(m,n,:)),squeeze(sc_NorCPM1(m,n,:)));
       R_sc2(m,n)=corr(squeeze(sc_EN4(m,n,:)),squeeze(sc_NorCPM2(m,n,:)));
       R_scf(m,n)=corr(squeeze(sc_EN4(m,n,:)),squeeze(sc_free(m,n,:)));
       R_sc3(m,n)=corr(squeeze(sc_EN4(m,n,:)),squeeze(sc_NorCPM1(m,n,:)+sc_NorCPM2(m,n,:))/2);
%       rmse_hc1(m,n)=sqrt(nanmean((hc_EN4(m,n,:)-hc_NorCPM1(m,n,:)).^2)/DEPTH_LIM^2-bias_hc1(m,n)^2);
%       rmse_hc2(m,n)=sqrt(nanmean((hc_EN4(m,n,:)-hc_NorCPM2(m,n,:)).^2)/DEPTH_LIM^2-bias_hc2(m,n)^2);
%       rmse_hcf(m,n)=sqrt(nanmean((hc_EN4(m,n,:)-hc_free(m,n,:)).^2)/DEPTH_LIM^2-bias_hcf(m,n)^2);
 
%       rmse_sc1(m,n)=sqrt(nanmean((sc_EN4(m,n,:)-sc_NorCPM1(m,n,:)).^2)/DEPTH_LIM^2-bias_sc1(m,n)^2);
%       rmse_sc2(m,n)=sqrt(nanmean((sc_EN4(m,n,:)-sc_NorCPM2(m,n,:)).^2)/DEPTH_LIM^2-bias_sc2(m,n)^2);
%       rmse_scf(m,n)=sqrt(nanmean((sc_EN4(m,n,:)-sc_free(m,n,:)).^2)/DEPTH_LIM^2-bias_scf(m,n)^2);
     end
   end
   save(['Correlations_hc_sc_yearly_' num2str(DEPTH_LIM) 'm_EN4_NorCPM.mat'],'R_hc1','R_hc2','R_hcf','R_hc3','R_sc1','R_sc2','R_scf','R_sc3','bias_hc1','bias_hc2','bias_hcf','bias_sc1','bias_sc2','bias_scf');
else
   load(['Correlations_hc_sc_yearly_' num2str(DEPTH_LIM) 'm_EN4_NorCPM.mat'])
end

close all
if 1
%%%%%%%%%%%%%%%%%%
figure(7)
clf
a=bias_hc1;
set(gcf, 'Renderer', 'opengl')
set(gcf, 'InvertHardCopy', 'off');
set(gca,'xticklabels',[],'yticklabels',[]);
%set(gca,'XAxisLocation','top','fontsize',16,'xtick',[]);
whitebg('w');
hold on
m_proj('Equidistant Cylindrical','lon',[-325 25],'lat',[-90 90]);
P1=m_pcolor(plon,plat,a)
set(P1,'LineStyle','none')
P1=m_pcolor(plon-360,plat,a);
set(P1,'LineStyle','none')
m_grid('xticklabels',[],'yticklabels',[]);
m_ungrid;
m_coast('patch',[0.5 0.5 0.5]);
colorbar('h','FontSize', 16)
colormap(anom(20));
caxis([-3 3]);
%title('SC corr between EN4 and reanalysis (SSTA)', 'fontsize', 20)
%title('SC ACC','fontsize', 20)
%print('-djpeg95',['Correlation_sc_anom_yearly_' num2str(DEPTH_LIM) '_EN4_free_fanf.jpg'])
%print('-depsc2',['Correlation_sc_anom_yearly_' num2str(DEPTH_LIM) '_EN4_fanf.eps'])
print('-depsc2',['Bias_hc_yearly_' num2str(DEPTH_LIM) '_i1.eps'])
%%%%%%%%%%%%%%%%%%
end

if 1
%%%%%%%%%%%%%%%%%%
figure(10)
clf
a=bias_sc1;
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 25],'lat',[-90 90]);
P1=m_pcolor(plon,plat,a)
set(P1,'LineStyle','none')
P1=m_pcolor(plon-360,plat,a);
set(P1,'LineStyle','none')
m_grid('xticklabels',[],'yticklabels',[]);
m_ungrid;
m_coast('patch',[0.5 0.5 0.5]);
colorbar('h','FontSize', 16)
colormap(anom(20));
caxis([-1 1])
%title('Reanalysis (SST, T and S) - reanalysis (SST)', 'FontSize', 20);
%print('-djpeg95',['Diff_Correlation_hc_detrend_anom_yearly_' num2str(DEPTH_LIM) '_NorCPM-free_fanf.jpg'])
%print('-depsc2',['Diff_Correlation_hc_detrend_anom_yearly_' num2str(DEPTH_LIM) '_yiguo-fanf.eps'])
print('-depsc2',['Bias_sc_yearly_' num2str(DEPTH_LIM) '_i1.eps'])
%%%%%%%%%%%%%%%%%%
end

if 1
figure(19)
clf
a=abs(bias_hc1)-abs(bias_hcf);
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 25],'lat',[-90 90]);
P1=m_pcolor(plon,plat,a);
set(P1,'LineStyle','none');
P1=m_pcolor(plon-360,plat,a);
set(P1,'LineStyle','none');
m_grid('xticklabels',[],'yticklabels',[]);
m_ungrid;
m_coast('patch',[0.5 0.5 0.5]);
colorbar('h','FontSize', 16)
colormap(anom(20));
caxis([-0.4 0.4]);
%title('SC corr between EN4 and reanalysis (SSTA)', 'fontsize', 20)
%title('SC ACC','fontsize', 20)
%print('-djpeg95',['Correlation_sc_anom_yearly_' num2str(DEPTH_LIM) '_EN4_free_fanf.jpg'])
%print('-depsc2',['Correlation_sc_anom_yearly_' num2str(DEPTH_LIM) '_EN4_fanf.eps'])
print('-depsc2',['Bias_hc_yearly_' num2str(DEPTH_LIM) '_i1-free.eps'])

figure(20)
clf
a=abs(bias_hc2)-abs(bias_hc1);
set(gcf, 'Renderer', 'opengl')
set(gca,'xticklabels',[],'yticklabels',[]);
set(gcf, 'InvertHardCopy', 'off');
%set(gca,'XAxisLocation','top','fontsize',16);
whitebg('w');
hold on
m_proj('Equidistant Cylindrical','lon',[-325 25],'lat',[-90 90]);
P1=m_pcolor(plon,plat,a);
set(P1,'LineStyle','none');
P1=m_pcolor(plon-360,plat,a);
set(P1,'LineStyle','none');
m_grid('xticklabels',[],'yticklabels',[]);
m_ungrid;
m_coast('patch',[0.5 0.5 0.5]);
colorbar('h','FontSize', 16);
colormap(anom(20));
caxis([-0.4 0.4]);
%title('SC corr between EN4 and reanalysis (SSTA)', 'fontsize', 20)
%title('SC ACC','fontsize', 20)
%print('-djpeg95',['Correlation_sc_anom_yearly_' num2str(DEPTH_LIM) '_EN4_free_fanf.jpg'])
%print('-depsc2',['Correlation_sc_anom_yearly_' num2str(DEPTH_LIM) '_EN4_fanf.eps'])
print('-depsc2',['Bias_hc_yearly_' num2str(DEPTH_LIM) '_i2-i1.eps'])

figure(21)
clf
a=abs(bias_sc1)-abs(bias_scf);
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 25],'lat',[-90 90]);
P1=m_pcolor(plon,plat,a);
set(P1,'LineStyle','none');
P1=m_pcolor(plon-360,plat,a);
set(P1,'LineStyle','none');
m_grid('xticklabels',[],'yticklabels',[]);
m_ungrid;
m_coast('patch',[0.5 0.5 0.5]);
colorbar('h','FontSize', 16);
colormap(anom(20));
caxis([-0.1 0.1]);
%title('SC corr between EN4 and reanalysis (SSTA)', 'fontsize', 20)
%title('SC ACC','fontsize', 20)
%print('-djpeg95',['Correlation_sc_anom_yearly_' num2str(DEPTH_LIM) '_EN4_free_fanf.jpg'])
%print('-depsc2',['Correlation_sc_anom_yearly_' num2str(DEPTH_LIM) '_EN4_fanf.eps'])
print('-depsc2',['Bias_sc_yearly_' num2str(DEPTH_LIM) '_i1-free.eps'])

figure(22)
clf
a=abs(bias_sc2)-abs(bias_sc1);
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 25],'lat',[-90 90]);
P1=m_pcolor(plon,plat,a);
set(P1,'LineStyle','none');
P1=m_pcolor(plon-360,plat,a);
set(P1,'LineStyle','none');
m_grid('xticklabels',[],'yticklabels',[]);
m_ungrid;
m_coast('patch',[0.5 0.5 0.5]);
colorbar('h','FontSize', 16);
colormap(anom(20));
caxis([-0.1 0.1]);
%title('SC corr between EN4 and reanalysis (SSTA)', 'fontsize', 20)
%title('SC ACC','fontsize', 20)
%print('-djpeg95',['Correlation_sc_anom_yearly_' num2str(DEPTH_LIM) '_EN4_free_fanf.jpg'])
%print('-depsc2',['Correlation_sc_anom_yearly_' num2str(DEPTH_LIM) '_EN4_fanf.eps'])
print('-depsc2',['Bias_sc_yearly_' num2str(DEPTH_LIM) '_i2-i1.eps'])
end