function [corr_norcpm,corr_free, sig_norcpm, sig_free]=Validation_SIC_decadalhindcast
%=== info on output: quote corr_est_df.m
% this functions calculates the correlation coefficient and the probability that the sample correlation is r, 
% even though the true correlation is 0. 
% The difference to the built-in matlab function corr is the calculation of the 
% degrees of freedom. Matlab uses df=n-2, here the effective degree of freedom is calculated
%
% Zur interpretation von p_val: wenn p_val < 0.05 dann lehnen wir Nullhypothese r=0 ab und die correlation ist
% signifikant unterschiedlich zu 0. STIMMT NICHT t-Test nach von Storch und Zwiers (p.149 eq. 8.7)
addpath('/cluster/home/madimm/matlab')
addpath('/cluster/home/madimm/matlab/m_map')
addpath('/cluster/home/madimm/matlab/Colormap')
%==time
yr_start=1960
yr_end  =2018
nb_yr=yr_end-yr_start + 1;

determine_free=1;
%==model data 
DPATH='/tos-project4/NS9039K/shared/norcpm/cases/NorESM1-CMIP6/';
FREE_PATH = [DPATH 'noresm1-cmip6_historical_19500101/Ensemble_mean/ice/']; %'/tos-project4/NS9039K/shared/norcpm/cases/NorESM1-ME_historicalExt_noAssim/';
if (0)
  CASE='i1';
  ASSIM_PATH = [DPATH 'norcpm-cmip6_hindcast/Processed/']; % '/tos-project4/NS9039K/shared/norcpm/cases/NorCPM/NorCPM_V1/ana_19800115_me_hindcasts/Processed/';
  CASENAME='norcpm-cmip6_hindcast_';
  OUTPUTPATH='/cluster/home/madimm/PROGRAMS_PaperIngoDCPP/PICS/';
else
  CASE='i2';
  ASSIM_PATH = [DPATH 'noresm1-cmip6_hindcast/Processed/']; %'/tos-project4/NS9039K/shared/norcpm/cases/NorCPM/NorCPM_V2a/hindcast/Processed/';
  OUTPUTPATH='/cluster/home/madimm/PROGRAMS_PaperIngoDCPP/PICS/';
  CASENAME='noresm1-cmip6_hindcast_';
end;
%=================
%== GRID
idm=320; jdm=384;
grid_path='/cluster/home/madimm/norcpm_link/SFE/';
plon=ncgetvar([grid_path 'grid.nc'],'plon');   plat=ncgetvar([grid_path 'grid.nc'],'plat');
parea=ncgetvar([grid_path 'grid.nc'],'parea'); mask=find(plat<0); parea(mask)=0;

% pivot points
ipiv=int16(ncread('/cluster/projects/nn9039k/NorCPM/Old/Input/NorESM/f19_g16/HADISST2/pivots_hadi-micom.nc','ipiv'));
jpiv=int16(ncread('/cluster/projects/nn9039k/NorCPM/Old/Input/NorESM/f19_g16/HADISST2/pivots_hadi-micom.nc','jpiv'));
%iipiv=int16(ncread('/cluster/projects/nn9039k/NorCPM/Old/Input/NorESM/f19_g16/HADISST2/pivots_SST.nc','ipiv'));
%jjpiv=int16(ncread('/cluster/projects/nn9039k/NorCPM/Old/Input/NorESM/f19_g16/HADISST2/pivots_SST.nc','jpiv'));
%==OBS + OBS grid
O_PATH='/cluster/projects/nn9039k/NorCPM/Obs/ICEC/HADISST2/HadISST_ice_1870_2018.nc'; %contains SIC data from 1870 Jan to 2018 Dec, HadISST1
dlon=360; dlat=180;

%===================================================================
%  given: monthly ensemble averages, create the time averages : yr1, yr2-3 and yr4-5
%
model=zeros(idm,jdm,nb_yr,3);
obser=zeros(idm,jdm,nb_yr,3);
obstmp=zeros(dlon,dlat,nb_yr,3);

if (determine_free)
  free=zeros(idm,jdm,nb_yr,3);
end; 


if (~exist([OUTPUTPATH CASE '_modelAICE_LYmeans.mat'],'file'))
  if (CASE=='i1')
    load /cluster/home/madimm/PROGRAMS_PaperIngoDCPP/PICS/i1_modelAICE_LYmeans.mat
  else
    load /cluster/home/madimm/PROGRAMS_PaperIngoDCPP/PICS/i2_modelAICE_LYmeans.mat
  end; 
else
%===  the prediction begins in october, thus lead year 1: to be averaged from november until november
 time=0;
 for year=yr_start:1:yr_end
   year
   time=time+1;
   %== annual mean LY1 ================
  % for mm=1:12
  %    yr1=year + floor((mm+10-1)./12); %actual (lead) month
  %    mm1=mod((mm+10-1),12)+1;
  %    %===
  %    filename=[ASSIM_PATH 'hindcast' num2str(year) '1015_ice_' num2str(yr1) '-' num2str(mm1,'%2.2d') '.nc'];
  %    tmp = ncgetvar(filename, 'aice'); %unit: %
  %    model(:,:,time,1)=model(:,:,time,1)+ tmp(:,:); 
  % end;
  % model(:,:,time,1)= model(:,:,time,1)./12.;
   %
   %== biannual mean LY2-3 ================
   for mm=13:24 %36 
      yr1=year + floor((mm+10-1)./12); %actual (lead) month
      mm1=mod((mm+10-1),12)+1;
      %===
      filename=[ASSIM_PATH 'hindcast' num2str(year) '1015_ice_' num2str(yr1) '-' num2str(mm1,'%2.2d') '.nc'];
      tmp = ncgetvar(filename, 'aice'); %unit: %
      model(:,:,time,2)=model(:,:,time,2)+ tmp(:,:); 
   end;
  % model(:,:,time,2)= model(:,:,time,2)./24.;
  model(:,:,time,2)= model(:,:,time,2)./12.;
  %
  %== biannual mean LY4-5 ================
  % for mm=37:48%72 
  %    yr1=year + floor((mm+10-1)./12); %actual (lead) month
  %    mm1=mod((mm+10-1),12)+1;
  %    %===
  %    filename=[ASSIM_PATH 'hindcast' num2str(year) '1015_ice_' num2str(yr1) '-' num2str(mm1,'%2.2d') '.nc'];
  %    tmp = ncgetvar(filename, 'aice'); %unit: %
  %    model(:,:,time,3)=model(:,:,time,3)+ tmp(:,:); 
  % end;
  % model(:,:,time,3)= model(:,:,time,3)./24.;
  %===
 end;
 save([OUTPUTPATH '/' CASE '_modelAICE_LYmeans_year2.mat'], 'model'); 
end;%if not exists


%========================= until here
%=== same with observations!
%1870 jan: first data
%1870 - 1959: 90 years (90*12)+1 = Jan 1960 
% first lead month is November, i.e. (90*12)+11
time=(90*12)+11;
obsice = ncgetvar(O_PATH,'sic');
obsice=obsice(:,:,time:end);
obsice=obsice.*100.0;
d=size(obsice);
%== annual mean LY1 ================ 
%time=0;
%for year=yr_start:1:yr_end
%   time=time+1;
%   nr=0.;
%   for mm=1:12
%    if (12*(time-1)+mm <=d(3))
%    obstmp(:,:,time,1)=obstmp(:,:,time,1)+ obsice(:,:,12*(time-1)+mm); 
%    nr=nr+1.;
%    end;
%   end;
%   obstmp(:,:,time,1)= obstmp(:,:,time,1)./nr;
%end;
%== biannual mean LY2 ================
time=0;
for year=yr_start:1:(yr_end-1)
   year
   time=time+1;
   nr=0.;
   for mm=13:48 
     if (12*(time-1)+mm <=d(3))
       obstmp(:,:,time,2)=obstmp(:,:,time,2)+ obsice(:,:,12*(time-1)+mm); 
       nr=nr+1;
     end;
   end;
   obstmp(:,:,time,2)= obstmp(:,:,time,2)./nr;
end;

%== biannual mean LY2-3 ================
%time=0;
%for year=yr_start:1:(yr_end-2)
%   year
%   time=time+1;
%   nr=0.;
%   for mm=13:36 
%     if (12*(time-1)+mm <=d(3))
%       obstmp(:,:,time,2)=obstmp(:,:,time,2)+ obsice(:,:,12*(time-1)+mm); 
%       nr=nr+1;
%     end;
%   end;
%   obstmp(:,:,time,2)= obstmp(:,:,time,2)./nr;
%end;
%== biannual mean LY4-5 ================
%time=0;
%for year=yr_start:1:(yr_end-5)
%   time=time+1;
%   nr = 0.;
%   for mm=37:72 
%    if (12*(time-1)+mm <=d(3))
%      nr=nr+1.;
%      obstmp(:,:,time,3)=obstmp(:,:,time,3)+ obsice(:,:,12*(time-1)+mm); 
%    end;
%   end;
%   obstmp(:,:,time,3)= obstmp(:,:,time,3)./nr;
%end;

%==project onto model grid
time=0;
for year=yr_start:1:yr_end
   time=time+1;
   for ii=1:3
      d2p=zeros(idm,jdm);
      nb =zeros(idm,jdm);
      for i=1:idm
 	for j=1:jdm
 	 d2p(i,j)=obstmp(ipiv(i,j),jpiv(i,j),time,ii);   
         nb(i,j)=nb(i,j)+1; 
        end
      end
      obser(:,:,time,ii)=obser(:,:,time,ii)+ d2p(:,:)./nb(:,:);
   end;
end; 
save([OUTPUTPATH '/obsAICE_LYmeans_LY2.mat'], 'obser'); 

if (determine_free)
  %time=0;
  %for year=yr_start:1:yr_end
     %year
     %time=time+1;
     %nr=0.;
     %== annual mean LY1 ================ 
     %for mm=1:12
     %   yr1=year + floor((mm+10-1)./12); %actual (lead) month
     %   mm1=mod((mm+10-1),12)+1;
     %   if (yr1<=2018)
     %     filename=[FREE_PATH 'free-average-ice' num2str(yr1) '-' num2str(mm1,'%2.2d') '.nc'];
     %     tmp = ncgetvar(filename, 'aice'); %unit: %
     %     free(:,:,time,1)=free(:,:,time,1)+ tmp(:,:); 
     %     nr=nr+1.;
     %   end;
     %end;
     %free(:,:,time,1)=free(:,:,time,1)./nr;
   %end;
   time=0;
   for year=yr_start:1:yr_end
     year
     nr=0.;time=time+1;
     %== biannual mean LY2-3 ================
     for mm=13:48%36 
        yr1=year + floor((mm+10-1)./12); %actual (lead) month
        mm1=mod((mm+10-1),12)+1;
        if (yr1<=2018)
          filename=[FREE_PATH 'free-average-ice' num2str(yr1) '-' num2str(mm1,'%2.2d') '.nc'];
          tmp = ncgetvar(filename, 'aice'); %unit: %
          free(:,:,time,2)=free(:,:,time,2)+ tmp(:,:); 
          nr=nr+1.;
        end;
     end;
     free(:,:,time,2)=free(:,:,time,2)./nr;
   end;
   %time=0;
   %for year=yr_start:1:yr_end
   %  year
   %  nr=0.;time=time+1;
   % %== biannual mean LY4-5 ================
   % for mm=37:72 
   %     yr1=year + floor((mm+10-1)./12); %actual (lead) month
   %     mm1=mod((mm+10-1),12)+1;
   %     if (yr1<=2018)
   %       filename=[FREE_PATH 'free-average-ice' num2str(yr1) '-' num2str(mm1,'%2.2d') '.nc'];
   %       tmp = ncgetvar(filename, 'aice'); %unit: %
   %       free(:,:,time,3)=free(:,:,time,3)+ tmp(:,:); 
   %       nr=nr+1.;
   %     end;
   %  end;
   %  free(:,:,time,3)=free(:,:,time,3)./nr;
   %end;
  save([OUTPUTPATH '/freeAICE_LYmeans_LY2.mat'], 'free'); 

end; 

%==== determine correlations and significance 
%     no bias correction
corr_norcpm=zeros(idm,jdm,3);
pval_norcpm=ones(idm,jdm,3);
ttest_norcpm=zeros(idm,jdm,3);
x=1:nb_yr;xx=reshape(x,[nb_yr,1]);
x2=1:nb_yr-2;xx2=reshape(x2,[nb_yr-2,1]);
x5=1:nb_yr-5;xx5=reshape(x5,[nb_yr-5,1]);

for i=1:idm
  for j=1:jdm
    % LY1
    am=reshape(model(i,j,:,1), [nb_yr,1]);
    ao=reshape(obser(i,j,:,1), [nb_yr,1]);
    %remove 2nd order regression
    regm=polyfit(xx, am, 2);
    rego=polyfit(xx, ao, 2);
    if (isnan(am(1))==0)
      [corr_norcpm(i,j,1), pval_norcpm(i,j,1)]=corr_est_df(ao-polyval(rego, xx),am-polyval(regm, xx));
      if (pval_norcpm(i,j,1)<=0.05)
       ttest_norcpm(i,j,1)=1;
      end; 
    end;
    % LY2-3
    am=reshape(model(i,j,1:end-2,2), [nb_yr-2,1]);
    ao=reshape(obser(i,j,1:end-2,2), [nb_yr-2,1]);
    %remove 2nd order regression
    regm=polyfit(xx2, am, 2);
    rego=polyfit(xx2, ao, 2);
    if (isnan(am(1))==0)
      [corr_norcpm(i,j,2), pval_norcpm(i,j,2)]=corr_est_df(ao-polyval(rego, xx2),am-polyval(regm, xx2));
      if (pval_norcpm(i,j,2)<=0.05)
       ttest_norcpm(i,j,2)=1;
      end; 
    end;
    % LY4-5
    am=reshape(model(i,j,1:end-5,3), [nb_yr-5,1]);
    ao=reshape(obser(i,j,1:end-5,3), [nb_yr-5,1]);
    %remove 2nd order regression
    regm=polyfit(xx5, am, 2);
    rego=polyfit(xx5, ao, 2);
    if (isnan(am(1))==0)
      [corr_norcpm(i,j,3), pval_norcpm(i,j,3)]=corr_est_df(ao-polyval(rego, xx5),am-polyval(regm, xx5));
      if (pval_norcpm(i,j,3)<=0.05)
       ttest_norcpm(i,j,3)=1;
      end; 
    end;
  end
end
%
 save([OUTPUTPATH '/' CASE '_CORR_AICE_LYmeans_LY2_2poly.mat'], 'corr_norcpm', 'pval_norcpm'); 

if (determine_free)
corr_free=zeros(idm,jdm,3);
pval_free=ones(idm,jdm,3);
ttest_free=zeros(idm,jdm,3);
x=1:nb_yr;xx=reshape(x,[nb_yr,1]);
x2=1:nb_yr-2;xx2=reshape(x2,[nb_yr-2,1]);
x5=1:nb_yr-5;xx5=reshape(x5,[nb_yr-5,1]);

for i=1:idm
  for j=1:jdm
    % LY1
    am=reshape(free(i,j,:,1), [nb_yr,1]);
    ao=reshape(obser(i,j,:,1), [nb_yr,1]);
   %remove 2nd order regression
    regm=polyfit(xx, am, 2);
    rego=polyfit(xx, ao, 2);
    if (isnan(am(1))==0)
      [corr_free(i,j,1), pval_free(i,j,1)]=corr_est_df(ao-polyval(rego, xx),am-polyval(regm, xx));
      if (pval_free(i,j,1)<=0.05)
       ttest_free(i,j,1)=1;
      end; 
    end;
    % LY2-3
    am=reshape(free(i,j,1:end-2,2), [nb_yr-2,1]);
    ao=reshape(obser(i,j,1:end-2,2), [nb_yr-2,1]);
   %remove 2nd order regression
    regm=polyfit(xx2, am, 2);
    rego=polyfit(xx2, ao, 2);
    if (isnan(am(1))==0)
      [corr_free(i,j,2), pval_free(i,j,2)]=corr_est_df(ao-polyval(rego, xx2),am-polyval(regm, xx2));
      if (pval_free(i,j,2)<=0.05)
       ttest_free(i,j,2)=1;
      end; 
    end;
    % LY4-5
    am=reshape(free(i,j,1:end-5,3), [nb_yr-5,1]);
    ao=reshape(obser(i,j,1:end-5,3), [nb_yr-5,1]);
    %remove 2nd order regression
    regm=polyfit(xx5, am, 2);
    rego=polyfit(xx5, ao, 2);
    if (isnan(am(1))==0)
      [corr_free(i,j,3), pval_free(i,j,3)]=corr_est_df(ao-polyval(rego, xx5),am-polyval(regm, xx5));
      if (pval_free(i,j,3)<=0.05)
       ttest_free(i,j,3)=1;
      end; 
    end;
  end
end;
save([OUTPUTPATH '/FREE_CORR_AICE_LYmeans_LY2_2poly.mat'], 'corr_free', 'pval_free'); 


end;
%======== 
% save correlations

%===================================================================
%                               plot    A) ARCTIC B) ANTARCTIC
%===================================================================
cd /cluster/home/madimm/norcpm_link/SFE/

%m_proj('stereographic','lon',-47,'lat',80.0,'rad',42,'rec','off') %Atlantic
%m_proj('stereographic','lon',-40,'lat',89.9,'rad',45,'rec','off');

m_proj('stereographic','lon',17,'lat',89.0,'rad',45,'rec','off');%Arctic
m_proj('stereographic','lon',17,'lat',-89.0,'rad',40,'rec','off');%Antarctic

%LY1
%=======
corre=corr_free(:,:,1);
ttitle= ['corr(SIC), @LY1 hindcast-' CASE ]; 
outputname=[OUTPUTPATH 'corrSIC_AntArctic_LY1_historical_poly2.eps'];

%outputname=[OUTPUTPATH 'corrSIC_Arctic_LY1_' CASE 'poly2.eps'];
%outputname=[OUTPUTPATH 'corrSIC_AntArctic_LY1_' CASE '.eps'];

corre=corr_norcpm(:,:,1);
ttitle= ['corr(SIC), @LY1 hindcast-' CASE ]; 
outputname=[OUTPUTPATH 'corrSIC_AntArctic_LY1_' CASE 'poly2.eps'];
%outputname=[OUTPUTPATH 'corrSIC_AntArctic_LY1_' CASE 'poly2.eps'];
close all; 
cnt=1;
fig=figure(cnt); cnt=cnt+1;
set(fig, 'units', 'inches', 'position', [0 0 10 10]);
micom_flat(corre,[.90 .90 .90]); m_grid;
% pclon1=ncgetvar('grid.nc','pclon');
% pclat1=ncgetvar('grid.nc','pclat');
%ice_index= aice_m(:,:,1)*0.0;
%ice_index((aice_m(:,:,1)==0) & (aice_o(:,:,1)==0))=1;
%pclon1=pclon1((ttest_norcpm==0) & (isnan(aice_m(:,:,1))==0) & (ice_index==0));
%pclat1=pclat1((ttest_MAM==0) & (isnan(aice_m(:,:,1))==0)& (ice_index==0));
%pclon1=pclon1((ttest_norcpm(:,:,1)==0));
%pclat1=pclat1((ttest_norcpm(:,:,1)==0));
%sss=size(pclon1)
%pclon2a=pclon1(1:5:floor(sss(1)/1),1:2:end );
%pclat2a=pclat1(1:5:floor(sss(1)/1),1:2:end );
%pclon2b=pclon1(pclat1<72 );
%pclat2b=pclat1(pclat1<72 );
%hold on;
%m_plot(pclon2a,pclat2a,'.k','markersize',8);
%m_plot(pclon2b(1:2:end,1:2:end),pclat2b(1:2:end,1:2:end),'.k','markersize',8);
%hold off;
title(ttitle, 'fontsize', 30); 
colormap(anomisuperwide(20)); colorbar('h'); h=colorbar('h','FontSize',26); 
caxis([-1 1]); m_coast('patch',[0.75 0.75 0.75]); set(gcf,'color','white')
print(gcf,'-depsc2',outputname);

%=========
%LY2-3
corre=corr_free(:,:,2);
%ttitle= ['corr(SIC), @LY2-3 hindcast-' CASE ];
ttitle= ['corr(SIC), @LY2-3 historical' ];  
%outputname=[OUTPUTPATH 'corrSIC_Arctic_LY23_' CASE '.eps'];
%outputname=[OUTPUTPATH 'corrSIC_AntArctic_LY23_' CASE '.eps'];
outputname=[OUTPUTPATH 'corrSIC_AntArctic_LY2_historical_poly2.eps'];

corre=corr_norcpm(:,:,2);
%ttitle= ['corr(SIC), @LY2-3 hindcast-' CASE ];
ttitle= ['corr(SIC), @LY2 hindcast-' CASE ];  
%outputname=[OUTPUTPATH 'corrSIC_Arctic_LY23_' CASE '.eps'];
outputname=[OUTPUTPATH 'corrSIC_AntArctic_LY2_' CASE 'poly2.eps'];
fig=figure(cnt); cnt=cnt+1;
set(fig, 'units', 'inches', 'position', [0 0 10 10]);
micom_flat(corre,[.90 .90 .90]); m_grid; title(ttitle, 'fontsize', 30);
colormap(anomisuperwide(20)); colorbar('h'); h=colorbar('h','FontSize',26); 
caxis([-1 1]); m_coast('patch',[0.75 0.75 0.75]); set(gcf,'color','white')
print(gcf,'-depsc2',outputname);


%LY4-5
corre=corr_free(:,:,3);
ttitle= ['corr(SIC), @LY4-5 ' CASE ]; 
%outputname=[OUTPUTPATH 'corrSIC_Arctic_LY45_' CASE '.eps'];
outputname=[OUTPUTPATH 'corrSIC_AntArctic_LY45_' CASE '.eps'];
outputname=[OUTPUTPATH 'corrSIC_AntArctic_LY45_historical_poly2.eps'];

corre=corr_norcpm(:,:,3);
ttitle= ['corr(SIC), @LY4-5 dcppA-assim-' CASE ]; 
%outputname=[OUTPUTPATH 'corrSIC_Arctic_LY45_' CASE '.eps'];
outputname=[OUTPUTPATH 'corrSIC_AntArctic_LY45_' CASE 'poly2.eps'];
fig=figure(cnt); cnt=cnt+1;
set(fig, 'units', 'inches', 'position', [0 0 10 10]);
micom_flat(corre,[.90 .90 .90]); m_grid; title(ttitle, 'fontsize', 30);
colormap(anomisuperwide(20)); colorbar('h'); h=colorbar('h','FontSize',26); 
caxis([-1 1]); m_coast('patch',[0.75 0.75 0.75]); set(gcf,'color','white')
print(gcf,'-depsc2',outputname);





end




%=================================
% OUTPUT ARCTIC and SOUTHERN OCEAN AVERAGE FOR LY1, LY2, 
% 1) FULLFIELD
%
% a) FREE
%
% b) persistence
%
% c) hindcast i1
%
% d) hindcast i2
%

% 2) 2nd order regression removed
%
% a) FREE
%
% b) persistence
%
% c) hindcast i1
%
% d) hindcast i2
%


%=================================
% === shell commands for averaging
%if (0)
%%creating ensemble averages
%module load  NCO/4.7.2-intel-2018a
%PATHNAME='/tos-project4/NS9039K/shared/norcpm/cases/NorESM1-CMIP6/'
%CASENAME='noresm1-cmip6_hindcast'
%cd ${PATHNAME}/${CASENAME}
%Y1=1960
%Y2=2018
%for year in `seq ${Y1} ${Y2}`
% do
% echo ${year}
%% for yr in `seq 00 10`
% do
%    mm=`expr $year + $yr`
%    for month in `seq -w 01 12`
%    do
%    FILE=./Processed/hindcast${year}1015_ice_${mm}-${month}.nc
%    if [ -f $FILE ]; then 
%      echo "$FILE exists"
%    else
%         echo "      $FILE"
%	 ncea ./${CASENAME}_${year}1015/${CASENAME}_${year}1015_mem??/ice/hist/${CASENAME}_${year}1015_mem??.cice.h.${mm}-${month}.nc tmp.nc
%	 mv tmp.nc ./Processed/hindcast${year}1015_ice_${mm}-${month}.nc
%    fi 
% done
% done
% done
%
%
%PATHNAME='/tos-project4/NS9039K/shared/norcpm/cases/NorESM1-CMIP6/'
%CASENAME='noresm1-cmip6_hindcast'
%cd ${PATHNAME}/${CASENAME}
%Y1=1960
%Y2=2018
%for year in `seq ${Y1} ${Y2}`
% do
% echo '====' ${CASENAME}_${year}1015
% ls -l ./Processed/hindcast${year}1015_ice* | wc -l
% done
%
%end;


