% initial development : Madlen Kimmritz (madlen.kimmritz@nersc.no), Nov 2019
% This script determines monthly SIA from northern 
% and southern hemisphere - separately for each member 
% of reanalysis/hindcast/observation(cont) product and stores the results in a 
% .mat-file for further use 
%
addpath('/cluster/home/madimm/matlab')
%==time for reanalysis
yr_start=1950; yr_end  =2018; nb_yr=yr_end-yr_start + 1;

%=================================================================================
%==reanalysis
DPATH='/tos-project4/NS9039K/shared/norcpm/cases/NorESM1-CMIP6/';
if (1)
  %norcpm-cmip6_analysis_19500115, norcpm-cmip6_hindcast (i1 == V1a analysis + hindcasts)
  CASE='i1';
  ASSIM_PATH = [DPATH 'norcpm-cmip6_analysis_19500115/']; % '/tos-project4/NS9039K/shared/norcpm/cases/NorCPM/NorCPM_V1/ana_19800115_me_hindcasts/Processed/';
  HIND_PATH = [DPATH 'norcpm-cmip6_hindcast/norcpm-cmip6_hindcast_']; 
  HINDNAME='norcpm-cmip6_hindcast_'
  CASENAME='norcpm-cmip6_analysis_19500115_mem';
  OUTPUTPATH='/cluster/home/madimm/PROGRAMS_PaperIngoDCPP/PICS/';
else
  % noresm1-cmip6_analysis_19500115, noresm1-cmip6_hindcast (i2 == V1c analysis + hindcasts)  
  CASE='i2';
  ASSIM_PATH = [DPATH 'noresm1-cmip6_analysis_19500115/']; %'/tos-project4/NS9039K/shared/norcpm/cases/NorCPM/NorCPM_V2a/hindcast/Processed/';
  HIND_PATH = [DPATH 'noresm1-cmip6_hindcast/noresm1-cmip6_hindcast_']; 
  HINDNAME='noresm1-cmip6_hindcast_'
 OUTPUTPATH='/cluster/home/madimm/PROGRAMS_PaperIngoDCPP/PICS/';
  CASENAME='noresm1-cmip6_analysis_19500115_mem';
end;

SIA_Arctic=zeros(30,nb_yr*12); SIA_SOcean=zeros(30,nb_yr*12);

%=== 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'); 

areaN=parea*1.;areaN(:,1:180)=0.; %use upper half
areaS=parea*1.;areaS(:,180:end)=0.;


%======= SIA (North and SOUTH) for reanalysis
if (0)
for mem=1:1:30
  mem
  time=0;
  for year=yr_start:1:yr_end
    year
    for mm=1:12
      time=time+1;
      filename=[ASSIM_PATH CASENAME num2str(mem,'%2.2d') '/ice/hist/' CASENAME num2str(mem,'%2.2d') '.cice.h.' num2str(year) '-' num2str(mm,'%2.2d') '.nc'];
      tmp=ncgetvar(filename, 'aice'); %unit: %
      tmp=tmp./100.;
      AN=nansum(nansum(areaN.*tmp));
      AS=nansum(nansum(areaS.*tmp));
      SIA_Arctic(mem,time)=AN(1,1);
      SIA_SOcean(mem,time)=AS(1,1);
    end; 
  end; 
end; 
save([OUTPUTPATH '/SIA_' CASE '_time.mat'], 'SIA_Arctic', 'SIA_SOcean'); 
end;



%=================================================================================
%==== FREE: /projects/NS9039K/shared/norcpm/cases/NorESM1-CMIP6/noresm1-cmip6_historical (free historical ensemble)
  
  FREE_PATH = [DPATH 'noresm1-cmip6_historical_19500101/'];
  FREE_PATH1 = [DPATH 'noresm1-cmip6_historical_20150101/'];
  OUTPUTPATH='/cluster/home/madimm/PROGRAMS_PaperIngoDCPP/PICS/';
  CASENAMEF='noresm1-cmip6_historical_19500101_mem';
  CASENAMEF1='noresm1-cmip6_historical_20150101_mem';
  nbf=2014-1950+1;
  SIA_Arctic_FREE=zeros(30,nb_yr*12);
  SIA_SOcean_FREE=zeros(30,nb_yr*12);
for mem=1:1:30
  mem
  time=0;
  for year=yr_start:1:2014
    year
    for mm=1:12
      time=time+1;
      filename=[FREE_PATH CASENAMEF num2str(mem,'%2.2d') '/ice/hist/' CASENAMEF num2str(mem,'%2.2d') '.cice.h.' num2str(year) '-' num2str(mm,'%2.2d') '.nc'];
      tmp=ncgetvar(filename, 'aice'); %unit: %
      tmp=tmp./100.;
      AN=nansum(nansum(areaN.*tmp));
      AS=nansum(nansum(areaS.*tmp));
      SIA_Arctic_FREE(mem,time)=AN(1,1);
      SIA_SOcean_FREE(mem,time)=AS(1,1);
    end; 
  end; 
  %==
  for year=2015:1:yr_end
    year
    for mm=1:12
      time=time+1;
      filename=[FREE_PATH1 CASENAMEF1 num2str(mem,'%2.2d') '/ice/hist/' CASENAMEF1 num2str(mem,'%2.2d') '.cice.h.' num2str(year) '-' num2str(mm,'%2.2d') '.nc'];
      tmp=ncgetvar(filename, 'aice'); %unit: m, sea ice volume per unit area
      tmp=tmp./100.;
      AN=nansum(nansum(areaN.*tmp));
      AS=nansum(nansum(areaS.*tmp));
      SIA_Arctic_FREE(mem,time)=AN(1,1);
      SIA_SOcean_FREE(mem,time)=AS(1,1);
    end; 
  end; 

end; 

save([OUTPUTPATH '/SIA_FREE_time.mat'], 'SIA_Arctic_FREE', 'SIA_SOcean_FREE'); 



%=================================================================================
%======= SIA (North and SOUTH) for observations
if(0)
%==OBS + OBS grid
det_obs = 1
if (det_obs)
  O_PATH='/cluster/projects/nn9039k/NorCPM/Obs/ICEC/HADISST2/HadISST_ice_1870_2018.nc'; %SIC (1870 Jan - 2018 Dec, HadISST1)
  dlon=360; dlat=180;
  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'));
  SIA_HadISST_Arctic=zeros(1,nb_yr*12);
  SIA_HadISST_SOcean=zeros(1,nb_yr*12);
  tmpo=zeros(dlon,dlat);

  %1870 jan: first data
  %1870 - 1949: 80 years (80*12)+1 = Jan 1950 
  time=(80*12)+1;
  obsice = ncgetvar(O_PATH,'sic');  obsice=obsice(:,:,time:end); d=size(obsice);

  time=0;
  for year=yr_start:1:yr_end
    year
    for mm=1:12
      time=time+1;
      tmpo=obsice(:,:,time);
      %== transfer to model grid
      tmp=zeros(idm,jdm);
      nb =zeros(idm,jdm);
      for i=1:idm
 	for j=1:jdm
 	 tmp(i,j)=tmp(i,j) + tmpo(ipiv(i,j),jpiv(i,j));   
         nb(i,j)=nb(i,j)+1; 
        end
      end
      tmp=tmp./nb;
      %calc SIA 
      AN=nansum(nansum(areaN.*tmp));
      AS=nansum(nansum(areaS.*tmp));
      SIA_HadISST_Arctic(time)=AN(1,1);
      SIA_HadISST_SOcean(time)=AS(1,1);
    end; 
  end; 
  save([OUTPUTPATH '/SIA_HadISST_time.mat'], 'SIA_HadISST_Arctic', 'SIA_HadISST_SOcean'); 
end;
end;
%=================================================================================
% == hindcast (10 members, 10 years hindcast length)
yr_start=1960; yr_end  =2018; nb_yr=yr_end-yr_start + 1;
SIA_Arctic_hind=zeros(10,nb_yr, 10*12); SIA_SOcean_hind=zeros(10,nb_yr, 10*12);

for mem=1:10
  nyear=0;
  for year=yr_start:1:yr_end 
    nyear=nyear+1;
    CASENAME=[HINDNAME num2str(year)  '1015_mem' num2str(mem,'%2.2d')]
   for mm=1:120
      yr1=year + floor(((mm+10-1))./12); %lead month and year
      mm1=mod((mm+10-1),12)+1;
      %disp([year, mm, mm1, yr1]);      
      filename=[HIND_PATH num2str(year)  '1015/' CASENAME  '/ice/hist/' CASENAME '.cice.h.' num2str(yr1) '-' num2str(mm1,'%2.2d') '.nc'];
      %disp(filename);
      tmp=ncgetvar(filename, 'aice');
      tmp=tmp./100.;
      AN=0.;AS=0.;
      AN=nansum(nansum(areaN.*tmp));
      AS=nansum(nansum(areaS.*tmp));
      SIA_Arctic_hind(mem,nyear,mm)=AN(1,1);
      SIA_SOcean_hind(mem,nyear,mm)=AS(1,1);
   end;
  end; 
end; 
save([OUTPUTPATH '/SIA_' CASE 'hind_time.mat'], 'SIA_Arctic_hind', 'SIA_SOcean_hind'); 



