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

%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'
for j=1:dlat
   lat(j)=-84+j;
end
for i=1:dlon
   lon(i)=i;
end
for j=2:dlat-1
end
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(1);clf;hold on
fld=bias_hc1;
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.15 0.8 0.85])
set(gcf, 'InvertHardCopy', 'off');
whitebg('w');
hold on
P1=m_pcolor(lon,lat,fld)
set(P1,'LineStyle','none')
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
set(gca,'DataAspectRatio',[1.2 1 1e+16])
m_coast('patch',[0.5 0.5 0.5]);
colormap(anom(20));
caxis([-3 3]);
h=colorbar('southoutside','FontSize', 20,'xtick',-3:3);
%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_hc_yearly_' num2str(DEPTH_LIM) '_i1.eps'])
%%%%%%%%%%%%%%%%%%
end

if 1
%%%%%%%%%%%%%%%%%%
figure(2);clf;hold on
fld=bias_sc1;
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
P1=m_pcolor(lon,lat,fld)
set(P1,'LineStyle','none')
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
set(gca,'DataAspectRatio',[1.2 1 1e+16])
m_coast('patch',[0.5 0.5 0.5]);
colormap(anom(20));
caxis([-1 1]);
h=colorbar('southoutside','FontSize', 20,'xtick',-1:0.5:1);
shrink=0.75;
pos=get(h,'pos');
set(h,'pos',[pos(1)+(1-shrink)/2*pos(3)  0.4*pos(2) shrink*pos(3) shrink*pos(4)]);
print('-depsc2',['Bias_sc_yearly_' num2str(DEPTH_LIM) '_i1.eps'])
%%%%%%%%%%%%%%%%%%
end

if 1
%%%%%%%%%%%%%%%%%%
figure(3);clf;hold on
fld=abs(bias_hc1)-abs(bias_hcf);
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
P1=m_pcolor(lon,lat,fld)
set(P1,'LineStyle','none')
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
set(gca,'DataAspectRatio',[1.2 1 1e+16])
m_coast('patch',[0.5 0.5 0.5]);
colormap(anom(20));
caxis([-0.4 0.4]);
h=colorbar('southoutside','FontSize', 20,'xtick',-0.4:0.2:0.4);
shrink=0.75;
pos=get(h,'pos');
set(h,'pos',[pos(1)+(1-shrink)/2*pos(3)  0.4*pos(2) shrink*pos(3) shrink*pos(4)]);
print('-depsc2',['Bias_hc_yearly_' num2str(DEPTH_LIM) '_i1-free.eps'])
%%%%%%%%%%%%%%%%%%
end

if 1
%%%%%%%%%%%%%%%%%%
figure(4);clf;hold on
fld=abs(bias_sc1)-abs(bias_scf);
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
P1=m_pcolor(lon,lat,fld)
set(P1,'LineStyle','none')
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
set(gca,'DataAspectRatio',[1.2 1 1e+16])
m_coast('patch',[0.5 0.5 0.5]);
colormap(anom(20));
caxis([-0.1 0.1]);
h=colorbar('southoutside','FontSize', 20,'xtick',-0.1:0.05:0.1);
shrink=0.75;
pos=get(h,'pos');
set(h,'pos',[pos(1)+(1-shrink)/2*pos(3)  0.4*pos(2) shrink*pos(3) shrink*pos(4)]);
print('-depsc2',['Bias_sc_yearly_' num2str(DEPTH_LIM) '_i1-free.eps'])
%%%%%%%%%%%%%%%%%%
end

if 1
%%%%%%%%%%%%%%%%%%
figure(5);clf;hold on
fld=abs(bias_hc2)-abs(bias_hc1);
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
P1=m_pcolor(lon,lat,fld)
set(P1,'LineStyle','none')
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
set(gca,'DataAspectRatio',[1.2 1 1e+16])
m_coast('patch',[0.5 0.5 0.5]);
colormap(anom(20));
caxis([-0.4 0.4]);
h=colorbar('southoutside','FontSize', 20,'xtick',-0.4:0.2:0.4);
shrink=0.75;
pos=get(h,'pos');
set(h,'pos',[pos(1)+(1-shrink)/2*pos(3)  0.4*pos(2) shrink*pos(3) shrink*pos(4)]);
print('-depsc2',['Bias_hc_yearly_' num2str(DEPTH_LIM) '_i2-i1.eps'])
%%%%%%%%%%%%%%%%%%
end

if 1
%%%%%%%%%%%%%%%%%%
figure(6);clf;hold on
fld=abs(bias_sc2)-abs(bias_sc1);
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
P1=m_pcolor(lon,lat,fld)
set(P1,'LineStyle','none')
m_grid('xticklabels',[],'yticklabels',[],'linestyle','none');
set(gca,'DataAspectRatio',[1.2 1 1e+16])
m_coast('patch',[0.5 0.5 0.5]);
colormap(anom(20));
caxis([-0.1 0.1]);
h=colorbar('southoutside','FontSize', 20,'xtick',-0.1:0.05:0.1);
shrink=0.75;
pos=get(h,'pos');
set(h,'pos',[pos(1)+(1-shrink)/2*pos(3)  0.4*pos(2) shrink*pos(3) shrink*pos(4)]);
print('-depsc2',['Bias_sc_yearly_' num2str(DEPTH_LIM) '_i2-i1.eps'])
%%%%%%%%%%%%%%%%%%
end




