% initial development : Madlen Kimmritz (madlen.kimmritz@nersc.no), Nov 2019
% This script determines monthly SIV 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 
%
% model data: siv=hi*parea, as hi= ice volume per unit grid cell area
%
addpath('~/matlab')
%==time for reanalysis
yr_start=1950; yr_end  =2018; nb_yr=yr_end-yr_start + 1;

%=================================================================================
%==reanalysis
DPATH='/projects/NS9039K/shared/norcpm/cases/NorESM1-CMIP6/';
if (0)
  %norcpm-cmip6_analysis_19500115, norcpm-cmip6_hindcast (i1 == V1a analysis + hindcasts)
  CASE='i1';
  ASSIM_PATH = [DPATH 'norcpm-cmip6_analysis_19500115/'];
  HIND_PATH = [DPATH 'norcpm-cmip6_hindcast/norcpm-cmip6_hindcast_']; 
  HINDNAME='norcpm-cmip6_hindcast_'
  CASENAME='norcpm-cmip6_analysis_19500115_mem';
  OUTPUTPATH='/scratch/ingo/PROGRAMS_PaperIngoDCPP/PICS/';
else
  % noresm1-cmip6_analysis_19500115, noresm1-cmip6_hindcast (i2 == V1c analysis + hindcasts)  
  CASE='i2';
  ASSIM_PATH = [DPATH 'noresm1-cmip6_analysis_19500115/'];
  HIND_PATH = [DPATH 'noresm1-cmip6_hindcast/noresm1-cmip6_hindcast_']; 
  HINDNAME='noresm1-cmip6_hindcast_'
 OUTPUTPATH='/scratch/ingo/PROGRAMS_PaperIngoDCPP/PICS/';
  CASENAME='noresm1-cmip6_analysis_19500115_mem';
end;

SIV_Arctic=zeros(30,nb_yr*12); SIV_SOcean=zeros(30,nb_yr*12);

%=== GRID
idm=320; jdm=384;
grid_path='/projects/NS9039K/shared/ingo/jupyter/Data/';
plon=ncread([grid_path 'grid.nc'],'plon');   plat=ncread([grid_path 'grid.nc'],'plat');
parea=ncread([grid_path 'grid.nc'],'parea'); 

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


%======= SIV (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=ncread(filename, 'hi'); %unit: m, sea ice volume per unit area
      tmp(~isnan(tmp))=min(10,tmp(~isnan(tmp)));
      AN=nansum(nansum(areaN.*tmp));
      AS=nansum(nansum(areaS.*tmp));
      SIV_Arctic(mem,time)=AN(1,1);
      SIV_SOcean(mem,time)=AS(1,1);
    end; 
  end; 
end; 
save([OUTPUTPATH '/SIV_' CASE '_time.mat'], 'SIV_Arctic', 'SIV_SOcean'); 
end;


%=================================================================================
%==== FREE: /projects/NS9039K/shared/norcpm/cases/NorESM1-CMIP6/noresm1-cmip6_historical (free historical ensemble)
if 0 

  FREE_PATH = [DPATH 'noresm1-cmip6_historical_19500101/'];
  FREE_PATH1 = [DPATH 'noresm1-cmip6_historical_20150101/'];
  OUTPUTPATH='/scratch/ingo/PROGRAMS_PaperIngoDCPP/PICS/';
  CASENAMEF='noresm1-cmip6_historical_19500101_mem';
  CASENAMEF1='noresm1-cmip6_historical_20150101_mem';
  nbf=2014-1950+1;
  SIV_Arctic_FREE=zeros(30,nb_yr*12);
  SIV_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=ncread(filename, 'hi'); %unit: m, sea ice volume per unit area
      tmp(~isnan(tmp))=min(10,tmp(~isnan(tmp)));
      AN=nansum(nansum(areaN.*tmp));
      AS=nansum(nansum(areaS.*tmp));
      SIV_Arctic_FREE(mem,time)=AN(1,1);
      SIV_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=ncread(filename, 'hi'); %unit: m, sea ice volume per unit area
      tmp(~isnan(tmp))=min(10,tmp(~isnan(tmp)));
      AN=nansum(nansum(areaN.*tmp));
      AS=nansum(nansum(areaS.*tmp));
      SIV_Arctic_FREE(mem,time)=AN(1,1);
      SIV_SOcean_FREE(mem,time)=AS(1,1);
    end; 
  end; 

end; 

save([OUTPUTPATH '/SIV_FREE_time.mat'], 'SIV_Arctic_FREE', 'SIV_SOcean_FREE'); 

end

%=================================================================================
%======= SIA (North and SOUTH) for observations
if(0)
%=== PIOMAS DATA for ARCTIC
%PIOMAS MNTHLY DATA (unit: 1000 km^3)
SIV_PIOMAS_Arctic=zeros(2018-1979+1,13);
SIV_PIOMAS_Arctic(1,1:12) = [27.704 30.171 32.045 32.951 32.295 29.787 23.661 18.409 16.911 17.854 20.124 23.201]; %1979
SIV_PIOMAS_Arctic(2,1:12) = [26.530 29.154 31.128 32.235 31.807 29.152 22.865 17.790 16.316 17.325 19.424 22.451];
SIV_PIOMAS_Arctic(3,1:12) = [25.317 27.770 29.818 30.744 30.034 26.817 20.173 14.635 12.808 13.965 16.145 19.115];
SIV_PIOMAS_Arctic(4,1:12) = [ 22.633 25.486 27.653 28.968 28.302 25.616 19.665 14.834 13.508 14.911 17.766 21.066];
SIV_PIOMAS_Arctic(5,1:12) = [ 24.457 27.279 29.385 30.386 30.170 27.758 21.854 16.738 15.200 16.427 18.954 21.962];
SIV_PIOMAS_Arctic(6,1:12) = [ 24.733 27.209 29.151 30.330 29.816 27.093 20.941 16.148 14.635 15.664 18.097 21.292];
SIV_PIOMAS_Arctic(7,1:12) = [ 24.639 27.264 29.453 30.864 30.517 27.488 20.768 16.012 14.583 15.838 18.289 21.285];
SIV_PIOMAS_Arctic(8,1:12) = [ 24.678 27.518 29.770 30.934 30.456 27.861 21.953 17.293 16.076 17.378 19.857 22.738];
SIV_PIOMAS_Arctic(9,1:12) = [ 26.025 28.840 30.663 31.790 31.502 28.814 22.274 16.708 15.353 16.703 19.255 22.299];
SIV_PIOMAS_Arctic(10,1:12) = [ 25.471 28.132 30.283 31.200 30.213 27.230 21.098 16.328 14.987 16.246 18.925 22.125];
SIV_PIOMAS_Arctic(11,1:12) = [ 25.363 27.837 29.424 30.112 29.621 27.122 21.075 16.198 14.770 15.941 18.654 21.797];
SIV_PIOMAS_Arctic(12,1:12) = [ 24.967 27.632 29.393 29.907 28.905 25.494 19.515 14.885 13.815 15.218 18.258 21.522];
SIV_PIOMAS_Arctic(13,1:12) = [ 24.689 27.477 29.665 30.742 30.005 26.792 20.014 15.005 13.593 15.017 17.688 20.894];
SIV_PIOMAS_Arctic(14,1:12) = [ 24.183 26.906 28.680 29.642 29.259 26.849 20.899 15.915 15.078 16.580 19.157 22.328];
SIV_PIOMAS_Arctic(15,1:12) = [ 25.300 27.720 29.462 30.427 29.575 25.805 18.509 13.479 12.441 13.964 16.983 20.398];
SIV_PIOMAS_Arctic(16,1:12) = [ 23.853 26.673 28.714 29.735 29.312 26.552 20.030 14.895 13.855 15.196 17.940 21.211];
SIV_PIOMAS_Arctic(17,1:12) = [ 24.199 26.431 27.953 28.447 27.311 23.810 17.186 12.438 11.236 12.128 14.837 18.178];
SIV_PIOMAS_Arctic(18,1:12) = [ 21.648 24.474 26.387 27.445 27.255 24.866 19.282 14.896 13.949 15.394 17.398 20.152];
SIV_PIOMAS_Arctic(19,1:12) = [ 23.497 26.315 28.317 29.370 28.741 25.716 19.242 14.381 13.228 14.199 16.615 19.916];
SIV_PIOMAS_Arctic(20,1:12) = [ 23.472 26.355 28.300 29.409 28.823 25.361 18.577 13.286 11.626 12.812 15.640 18.951];
SIV_PIOMAS_Arctic(21,1:12) = [ 22.428 25.324 27.299 28.450 27.924 24.881 18.517 13.076 11.045 12.433 15.377 18.394];
SIV_PIOMAS_Arctic(22,1:12) = [ 21.712 24.297 26.174 27.153 26.689 23.881 17.429 12.422 11.084 12.313 14.901 17.952];
SIV_PIOMAS_Arctic(23,1:12) = [ 21.259 24.128 26.414 27.632 26.870 23.828 17.901 13.397 12.275 13.321 15.764 18.605];
SIV_PIOMAS_Arctic(24,1:12) = [ 21.837 24.793 26.655 27.428 26.822 23.685 17.164 12.127 10.846 11.890 14.587 17.844];
SIV_PIOMAS_Arctic(25,1:12) = [ 21.236 24.127 26.308 27.244 26.273 22.911 16.486 11.496 10.284 11.204 13.686 16.848];
SIV_PIOMAS_Arctic(26,1:12) = [ 20.050 22.655 24.809 25.750 25.308 22.706 16.577 11.515 10.037 11.247 13.966 17.238];
SIV_PIOMAS_Arctic(27,1:12) = [ 20.254 22.667 24.820 26.045 25.354 21.617 15.225 10.708  9.283 10.180 12.847 15.983];
SIV_PIOMAS_Arctic(28,1:12) = [ 19.307 22.036 24.051 25.105 24.337 20.904 14.692 10.404  9.113  9.810 12.263 15.039];
SIV_PIOMAS_Arctic(29,1:12) = [ 18.322 20.813 22.995 23.752 23.094 19.190 12.119  7.615  6.526  7.107 10.404 14.167];
SIV_PIOMAS_Arctic(30,1:12) = [ 18.515 21.527 23.823 24.983 24.122 20.576 14.154  9.201  7.250  8.267 11.700 15.134];
SIV_PIOMAS_Arctic(31,1:12) = [ 18.794 21.682 23.794 24.944 23.880 19.732 12.833  8.272  6.925  7.630 10.759 14.184];
SIV_PIOMAS_Arctic(32,1:12) = [ 17.667 20.581 23.082 24.097 22.223 17.142 10.244  5.927  4.742  6.198  9.480 12.925];
SIV_PIOMAS_Arctic(33,1:12) = [ 16.206 19.320 21.390 22.500 21.139 16.500  9.546  5.489  4.480  5.719  9.250 12.989];
SIV_PIOMAS_Arctic(34,1:12) = [ 16.886 19.585 21.924 23.119 21.714 16.002  9.264  4.964  3.787  5.001  8.217 12.133];
SIV_PIOMAS_Arctic(35,1:12) = [ 15.799 19.318 21.964 23.122 21.865 17.538 10.538  6.394  5.479  6.953 10.076 13.787];
SIV_PIOMAS_Arctic(36,1:12) = [ 17.406 19.850 21.802 22.935 21.913 17.680 11.951  8.150  6.972  8.159 11.481 15.074];
SIV_PIOMAS_Arctic(37,1:12) = [ 18.445 21.462 23.213 24.230 23.027 18.557 11.649  7.086  5.852  7.013 10.298 14.001];
SIV_PIOMAS_Arctic(38,1:12) = [ 17.185 19.592 21.524 22.459 21.026 16.493 10.257  5.941  4.530  5.511  7.834 11.206];
SIV_PIOMAS_Arctic(39,1:12) = [ 14.641 17.366 19.576 20.664 19.834 15.400  9.141  5.450  4.679  6.108  9.179 12.638];
SIV_PIOMAS_Arctic(40,1:12) = [ 16.022 18.639 20.844 22.245 21.429 17.206 10.295  6.203  5.081  5.776  9.426 13.172]; %2018
for i=1:40
SIV_PIOMAS_Arctic(i,13) =mean(SIV_PIOMAS_Arctic(i,1:12));
end;
SIV_PIOMAS_Arctic= SIV_PIOMAS_Arctic(:,13);
save([OUTPUTPATH '/SIV_PIOMAS_time.mat'], 'SIV_PIOMAS_Arctic'); 
end
%==========================================================
%=== combined data from cryosat and hadisst1 
%==OBS + OBS grid
det_obs = 0
if (det_obs)
  yr_starto=2002;
  yr_endo=2017;
  nobs=2017-2002+1;
  %=========== SIT DATA
  %SIT availability: (2002-2017), cryosat2 and envisat
  OSIT_PATH='/cluster/shared/noresm/norcpm/Obs/ICET/Merged/'; 
  ipivSIT=int16(ncread('/cluster/shared/noresm/norcpm/Input/NorESM/f19_g16/ICET/pivots_ICET.nc','ipiv'));
  jpivSIT=int16(ncread('/cluster/shared/noresm/norcpm/Input/NorESM/f19_g16/ICET/pivots_ICET.nc','jpiv'));

  dtlon=432; dtlat=432; tmpto=zeros(dtlon,dtlat);
  SIT_Arctic=zeros(dtlon,dtlat,nobs*12);
  SIC_Arctic=zeros(dtlon,dtlat,nobs*12);
  SIT_found=zeros(1,nobs*12); %1 when data available, 0 else, availability in winter months, we take average only over 
                              %those months in year, where data available (= good estimate of annual avg)
  time=0;
  for year=yr_starto:1:yr_endo
    year
    for mm=1:12
      time=time+1; 
      if exist([OSIT_PATH  num2str(year) '_' num2str(mm,'%2.2d') '.nc'], 'file') == 2
         SIT_found(1,time)=1; time
         tmpto=ncgetvar([OSIT_PATH  num2str(year) '_' num2str(mm,'%2.2d') '.nc'], 'sea_ice_thickness');
         tmpco=ncgetvar([OSIT_PATH  num2str(year) '_' num2str(mm,'%2.2d') '.nc'], 'sea_ice_concentration');
         tmpco=tmpco./100.;
         SIC_Arctic(:,:,time)= tmpco;
         SIT_Arctic(:,:,time)= tmpto;
        %
      end;
    end; 
  end; 


  %========= SIC DATA
  %SIC availability: (1870-2018), HadISST1
  OSIC_PATH='/cluster/projects/nn9039k/NorCPM/Obs/ICEC/HADISST2/HadISST_ice_1870_2018.nc'; 
  ipivSIC=int16(ncread('/cluster/projects/nn9039k/NorCPM/Old/Input/NorESM/f19_g16/HADISST2/pivots_hadi-micom.nc','ipiv'));
  jpivSIC=int16(ncread('/cluster/projects/nn9039k/NorCPM/Old/Input/NorESM/f19_g16/HADISST2/pivots_hadi-micom.nc','jpiv'));
  SIC_HadISST_Arctic=zeros(idm,jdm,nobs*12);   SIC_HadISST_Arctic_T=zeros(dtlon,dtlat,nobs*12); area_T=zeros(dtlon,dtlat);
  dlon=360; dlat=180; tmpo=zeros(dlon,dlat);

  %1870 - 2001: 132 years (132*12)+1 = Jan 2002 
  time1=(132*12)+1; obsice = ncgetvar(OSIC_PATH,'sic');  obsice=obsice(:,:,time1:end); d=size(obsice);
  time=0;
  for year=yr_starto:1:yr_endo
    year
    for mm=1:12
      time=time+1; tmpo=obsice(:,:,time);
      %== transfer to model grid
      nb =zeros(idm,jdm);
      for i=1:idm
 	for j=1:jdm
 	  SIC_HadISST_Arctic(i,j,time)=SIC_HadISST_Arctic(i,j,time) + tmpo(ipivSIC(i,j),jpivSIC(i,j));   
         nb(i,j)=nb(i,j)+1; 
        end
      end
       SIC_HadISST_Arctic(:,:,time)= SIC_HadISST_Arctic(:,:,time)./nb;
      %====
      nbt =zeros(dtlon,dtlat);
      for i=1:dtlon
 	for j=1:dtlat
 	  SIC_HadISST_Arctic_T(i,j,time)=SIC_HadISST_Arctic_T(i,j,time) + SIC_HadISST_Arctic(ipivSIT(i,j),jpivSIT(i,j),time);   
         nbt(i,j)=nbt(i,j)+1; 
        end
      end
       SIC_HadISST_Arctic_T(:,:,time)= SIC_HadISST_Arctic_T(:,:,time)./nbt;
    end; 
  end; 
  %===
  nbt =zeros(dtlon,dtlat);
  for i=1:dtlon
    for j=1:dtlat
      area_T(i,j)=area_T(i,j) + areaN(ipivSIT(i,j),jpivSIT(i,j));   
      nbt(i,j)=nbt(i,j)+1; 
    end
  end
  area_T(:,:)= area_T(:,:)./nbt;



  %======== SIV DATA
  time=0;
  SIV_obsmerged_Arctic=zeros(dtlon,dtlat,nobs*12);
  for year=yr_starto:1:yr_endo
    year
    for mm=1:12
      time=time+1; 
%      SIV_obsmerged_Arctic(:,:,time)=SIC_HadISST_Arctic_T(:,:,time).*SIT_Arctic(:,:,time);
      SIV_obsmerged_Arctic(:,:,time)=SIC_Arctic(:,:,time).*SIT_Arctic(:,:,time);
     % SIV_obsmerged_Arctic(:,:,time)=(SIT_Arctic(:,:,time));
      SIV_obsmerged_Arctic(:,:,time)=SIV_obsmerged_Arctic(:,:,time).*area_T;
    end; 
  end; 
  SIV_obsmerged_Arctic_time=zeros(1,nobs*12);
  time=0;
  for year=yr_starto:1:yr_endo
    year
    for mm=1:12
      time=time+1; 
      SIV_obsmerged_Arctic_time(1,time)=nansum(nansum(SIV_obsmerged_Arctic(:,:,time)));
    end; 
  end; 

SIV_obsmerged_Arctic_time_merged = SIV_obsmerged_Arctic_time.*1.;
SIV_obsmerged_Arctic_time_HadISST = SIV_obsmerged_Arctic_time.*1.;
  save([OUTPUTPATH '/SIV_OBSmerged_time.mat'], 'SIV_obsmerged_Arctic_time_merged', 'SIV_obsmerged_Arctic_time_HadISST', 'SIT_found'); 

end;
%end;
%=================================================================================
% == hindcast (10 members, 10 years hindcast length)
yr_start=1960; yr_end  =2018; nb_yr=yr_end-yr_start + 1;
SIV_Arctic_hind=zeros(10,nb_yr, 10*12); SIV_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=ncread(filename, 'hi');
      tmp(~isnan(tmp))=min(10,tmp(~isnan(tmp)));
      AN=0.;AS=0.;
      AN=nansum(nansum(areaN.*tmp));
      AS=nansum(nansum(areaS.*tmp));
      SIV_Arctic_hind(mem,nyear,mm)=AN(1,1);
      SIV_SOcean_hind(mem,nyear,mm)=AS(1,1);
   end;
  end; 
end; 
save([OUTPUTPATH '/SIV_' CASE 'hind_time.mat'], 'SIV_Arctic_hind', 'SIV_SOcean_hind'); 



