% initial development : Madlen Kimmritz (madlen.kimmritz@nersc.no), Nov 2019
% Plot time development of SIA (annual avgs) for
%  - reanalysis (avg, max, min)
%  - observation
%  - hindcast (avg,max,min) branched every n years (n=10)
%  - optional: drift adjustment
% for southern and northern hemisphere
%====
addpath('/cluster/home/madimm/matlab')
%==time for reanalysis
yr_start=1950; yr_end  =2018; nyrs=yr_end-yr_start + 1;

%=== read data
OUTPUTPATH='/cluster/home/madimm/PROGRAMS_PaperIngoDCPP/PICS/';
CASE ='i2';
%reanalysis: 'SIA_Arctic', 'SIA_SOcean', dimensions: (30,nb_yr*12);
load([OUTPUTPATH 'SIA_' CASE '_time.mat']);
%hindcast: 'SIA_Arctic_hind', 'SIA_SOcean_hind', dimensions: (10,nb_yr, 10*12)
load([OUTPUTPATH 'SIA_' CASE 'hind_time.mat']);
% observations 'SIA_HadISST_Arctic', 'SIA_HadISST_SOcean', dimensions: (1,nb_yr*12);
load([OUTPUTPATH 'SIA_HadISST_time.mat'])
%free: 'SIA_Arctic_FREE', 'SIA_SOcean_FREE', 
load([OUTPUTPATH 'SIA_FREE_time.mat']);

%=== annual avgs
%observations and reanalysis start from 1950 jan
%hindcast start from 1960, nov
nmths=828;
nyrs=69;
SIA_SO_ann=zeros(30,69);SIA_AR_ann=zeros(30,69);
SIA_SO_obs_ann=zeros(1,69);SIA_AR_obs_ann=zeros(1,69);
SIA_SO_ann_hind=zeros(10,59,9);SIA_AR_ann_hind=zeros(10,59,9);%start 10 years later in 1960
SIA_SO_ann_FREE=zeros(30,69);SIA_AR_ann_FREE=zeros(30,69);

for i=1:nyrs
  for mem=1:30
    SIA_SO_ann(mem,i)=mean(SIA_SOcean(mem,(i-1)*12+1:i*12));
    SIA_AR_ann(mem,i)=mean(SIA_Arctic(mem,(i-1)*12+1:i*12));
  end; 
  SIA_SO_obs_ann(1,i)=mean(SIA_HadISST_SOcean(1,(i-1)*12+1:i*12)); 
  SIA_AR_obs_ann(1,i)=mean(SIA_HadISST_Arctic(1,(i-1)*12+1:i*12)); 
end;


for i=1:nyrs
  for mem=1:30
    SIA_SO_ann_FREE(mem,i)=mean(SIA_SOcean_FREE(mem,(i-1)*12+1:i*12));
    SIA_AR_ann_FREE(mem,i)=mean(SIA_Arctic_FREE(mem,(i-1)*12+1:i*12));
  end; 
end;



for i=1:nyrs-10 
  %hindcast starts 10 years later !=====
  for mem=1:10
   for lyr=1:9
    start=(lyr-1)*12+3; 
    ende=start+11;
    SIA_SO_ann_hind(mem,i,lyr)=mean(SIA_SOcean_hind(mem,i,start:ende));
    SIA_AR_ann_hind(mem,i,lyr)=mean(SIA_Arctic_hind(mem,i,start:ende));
   end;
  end; 
end;

%=== determine min, max, avg
SIA_SO_min=zeros(1,69);SIA_AR_min=zeros(1,69);
SIA_SO_max=zeros(1,69);SIA_AR_max=zeros(1,69);
SIA_SO_avg=zeros(1,69);SIA_AR_avg=zeros(1,69);

SIA_SO_min_hind=zeros(59,9);SIA_AR_min_hind=zeros(59,9);
SIA_SO_max_hind=zeros(59,9);SIA_AR_max_hind=zeros(59,9);
SIA_SO_avg_hind=zeros(59,9);SIA_AR_avg_hind=zeros(59,9);


SIA_SO_min_FREE=zeros(1,69);SIA_AR_min_FREE=zeros(1,69);
SIA_SO_max_FREE=zeros(1,69);SIA_AR_max_FREE=zeros(1,69);
SIA_SO_avg_FREE=zeros(1,69);SIA_AR_avg_FREE=zeros(1,69);



for i=1:nyrs
    SIA_SO_min(1,i)=min(SIA_SO_ann(:,i));
    SIA_AR_min(1,i)=min(SIA_AR_ann(:,i));
    SIA_SO_max(1,i)=max(SIA_SO_ann(:,i));
    SIA_AR_max(1,i)=max(SIA_AR_ann(:,i));
    SIA_SO_avg(1,i)=mean(SIA_SO_ann(:,i));
    SIA_AR_avg(1,i)=mean(SIA_AR_ann(:,i));
end;

for i=1:nyrs
    SIA_SO_min_FREE(1,i)=min(SIA_SO_ann_FREE(:,i));
    SIA_AR_min_FREE(1,i)=min(SIA_AR_ann_FREE(:,i));
    SIA_SO_max_FREE(1,i)=max(SIA_SO_ann_FREE(:,i));
    SIA_AR_max_FREE(1,i)=max(SIA_AR_ann_FREE(:,i));
    SIA_SO_avg_FREE(1,i)=mean(SIA_SO_ann_FREE(:,i));
    SIA_AR_avg_FREE(1,i)=mean(SIA_AR_ann_FREE(:,i));
end;

for i=1:nyrs-10
  for ly=1:9
    SIA_SO_min_hind(i,ly)=min(SIA_SO_ann_hind(:,i,ly));
    SIA_AR_min_hind(i,ly)=min(SIA_AR_ann_hind(:,i,ly));
    SIA_SO_max_hind(i,ly)=max(SIA_SO_ann_hind(:,i,ly));
    SIA_AR_max_hind(i,ly)=max(SIA_AR_ann_hind(:,i,ly));
    SIA_SO_avg_hind(i,ly)=mean(SIA_SO_ann_hind(:,i,ly));
    SIA_AR_avg_hind(i,ly)=mean(SIA_AR_ann_hind(:,i,ly));
  end;
end;


%=======
% pick every 10 years and plot min,max,avg hull
records=[1,11,21,31,41,51]+1; %1960=11th entry,1970,1980,1990,2000,2010 + 1
targx=zeros(6,9);
for iii=1:6
  targx(iii,:)=records(iii)+10:records(iii)+18;
end; 
cnt=31; 
%=== plot
xaxis=linspace(yr_start,yr_end,yr_end-yr_start+1);
x=1: yr_end-yr_start+1;
%ARCTIC
figure(cnt); cnt = cnt+1; set(gcf,'Position',[100 100 2000 500]); hold on; 
hfree= plot(SIA_AR_avg_FREE./1e12, 'Color', [0.99 0.6 0.1] , 'LineWidth', 5);

x2 = [x, fliplr(x)]; inBetween = [SIA_AR_min./1e12, fliplr(SIA_AR_max./1e12)]; hhh=fill(x2, inBetween, [0.8 0.8 0.8]);
set(hhh,'EdgeColor','none'); set(hhh,'facealpha',.6);
hreana= plot(SIA_AR_avg./1e12, 'k' , 'LineWidth', 8); hobs=plot(SIA_AR_obs_ann./1e12, 'g' , 'LineWidth', 5, 'LineStyle', ':');

for iii=1:5
temp1=(SIA_AR_min_hind(records(1,iii)-1,:)./1e12); temp2=(SIA_AR_max_hind(records(1,iii)-1,:)./1e12);
x3 = [targx(iii,:), fliplr(targx(iii,:))]; inBetween = [temp1, fliplr(temp2)]; hh=fill(x3, inBetween, 'b'); 
set(hh,'facealpha',.2); set(hh,'EdgeColor','none');
end; 
iii=6
temp1=(SIA_AR_min_hind(records(1,iii)-1,1:8)./1e12); temp2=(SIA_AR_max_hind(records(1,iii)-1,1:8)./1e12);
x3 = [targx(iii,1:8), fliplr(targx(iii,1:8))]; inBetween = [temp1, fliplr(temp2)]; hh=fill(x3, inBetween, 'b'); 
set(hh,'facealpha',.2); set(hh,'EdgeColor','none');

hhind=plot(targx(1,:), squeeze(SIA_AR_avg_hind(records(1,1)-1,:)./1e12), 'b' , 'LineWidth', 8);
for iii=2:5
plot(targx(iii,:), squeeze(SIA_AR_avg_hind(records(1,iii)-1,:)./1e12), 'b' , 'LineWidth', 8);
end;
iii=6; plot(targx(iii,1:8), squeeze(SIA_AR_avg_hind(records(1,iii)-1,1:8)./1e12), 'b' , 'LineWidth', 8);

grid on;
set(gca,'XTick', [0 10 20 30 40 50 60]+1); set(gca,'XTickLabel', [1950 1960 1970 1980 1990 2000 2010]);
set(gca,'fontsize',20); ttitle= ['Arctic SIA ' CASE ]; title(ttitle, 'FontSize',24);
xlabel('years'); ylabel('SIA [mio km^2]'); xlim([1 69]); ylim([9 13]);
hleglines = [hfree(1) hreana(1) hhind(1) hobs(1)];
hleg = legend(hleglines,'historical', ['assim-' CASE],['hindcast-' CASE],'HadISST1');%,'Location','eastoutside', 'FontSize',18);
hleg.Location='southwest';
outputname=[OUTPUTPATH 'SIA_Arctic_time' CASE 'F1.eps'];
print(gcf,'-depsc2',outputname);


%SOUTHERN OCEAN
figure(cnt); cnt = cnt+1; 
set(gcf,'Position',[100 100 2000 500]); hold on;
hfree= plot(SIA_SO_avg_FREE./1e12, 'Color', [0.99 0.6 0.1] , 'LineWidth', 5);
x2 = [x, fliplr(x)]; inBetween = [SIA_SO_min./1e12, fliplr(SIA_SO_max./1e12)]; hhh=fill(x2, inBetween, [0.8 0.8 0.8]);
set(hhh,'EdgeColor','none'); set(hhh,'facealpha',.6);
hreana= plot(SIA_SO_avg./1e12, 'k' , 'LineWidth', 8); hobs=plot(SIA_SO_obs_ann./1e12, 'g' , 'LineWidth', 5, 'Linestyle', ':');
%

for iii=1:5
temp1=(SIA_SO_min_hind(records(1,iii)-1,:)./1e12); temp2=(SIA_SO_max_hind(records(1,iii)-1,:)./1e12);
x3 = [targx(iii,:), fliplr(targx(iii,:))]; inBetween = [temp1, fliplr(temp2)]; hh=fill(x3, inBetween, 'b'); 
set(hh,'facealpha',.2); set(hh,'EdgeColor','none'); 
end; 
iii=6
temp1=(SIA_SO_min_hind(records(1,iii)-1,1:8)./1e12); temp2=(SIA_SO_max_hind(records(1,iii)-1,1:8)./1e12);
x3 = [targx(iii,1:8), fliplr(targx(iii,1:8))]; inBetween = [temp1, fliplr(temp2)]; hh=fill(x3, inBetween, 'b'); 
set(hh,'facealpha',.2); 
set(hh,'EdgeColor','none')


hhind=plot(targx(1,:), squeeze(SIA_SO_avg_hind(records(1,1)-1,:)./1e12), 'b' , 'LineWidth', 8);
for iii=2:5
plot(targx(iii,:), squeeze(SIA_SO_avg_hind(records(1,iii)-1,:)./1e12), 'b' , 'LineWidth', 8);
end;
iii=6; plot(targx(iii,1:8), squeeze(SIA_SO_avg_hind(records(1,iii)-1,1:8)./1e12), 'b' , 'LineWidth', 8);

grid on;
set(gca,'XTick', [0 10 20 30 40 50 60]+1); set(gca,'XTickLabel', [1950 1960 1970 1980 1990 2000 2010]);
set(gca,'fontsize',20); ttitle= ['Southern Ocean SIA ' CASE ]; title(ttitle, 'FontSize',24);
xlabel('years'); ylabel('SIA [mio km^2]'); xlim([1 69]);
hleglines = [hfree(1) hreana(1) hhind(1) hobs(1)];
hleg = legend(hleglines,'historical', ['assim-' CASE], ['hindcast-' CASE],'HadISST1');
hleg.Location='northeast';
outputname=[OUTPUTPATH 'SIA_SouthOcean_time' CASE 'F1.eps'];
print(gcf,'-depsc2',outputname);



%=============================

driftadj=1;

if (driftadj)
%=== following YEager et al.: transform raw data into anomalies to climatological forecast, Y'_jt = Yjt - bar(Y_t)
% bar(Y_t): average of start years for given lead t       (= clim_t)
% Yjt: ens avg foecast from start year j at lead time t   (= SIA_avg_hind)
clim_t_AR=zeros(1,9);
clim_t_SO=zeros(1,9);
for ly=1:9
  clim_t_AR(1,ly) = mean(SIA_AR_avg_hind(:,ly));
  clim_t_SO(1,ly) = mean(SIA_SO_avg_hind(:,ly));
end;
SIA_SO_avg_hind_anom=zeros(59,9); 
SIA_AR_avg_hind_anom=zeros(59,9);
SIA_SO_min_hind_anom=zeros(59,9); 
SIA_AR_min_hind_anom=zeros(59,9);
SIA_SO_max_hind_anom=zeros(59,9); 
SIA_AR_max_hind_anom=zeros(59,9);
for i=1:nyrs-10
for ly=1:9
  SIA_SO_avg_hind_anom(i,ly)=SIA_SO_avg_hind(i,ly)-clim_t_SO(1,ly);
  SIA_AR_avg_hind_anom(i,ly)=SIA_AR_avg_hind(i,ly)-clim_t_AR(1,ly);
  SIA_SO_min_hind_anom(i,ly)=SIA_SO_min_hind(i,ly)-clim_t_SO(1,ly);
  SIA_AR_min_hind_anom(i,ly)=SIA_AR_min_hind(i,ly)-clim_t_AR(1,ly);
  SIA_SO_max_hind_anom(i,ly)=SIA_SO_max_hind(i,ly)-clim_t_SO(1,ly);
  SIA_AR_max_hind_anom(i,ly)=SIA_AR_max_hind(i,ly)-clim_t_AR(1,ly);
end; end; 
mean_AR=mean(SIA_AR_avg(1,11:end)); %over same period as hindcast
mean_SO=mean(SIA_SO_avg(1,11:end)); %over same period as hindcast
mean_ARF=mean(SIA_AR_avg_FREE(1,11:end)); %over same period as hindcast
mean_SOF=mean(SIA_SO_avg_FREE(1,11:end)); %over same period as hindcast

%===========
%=== plot
%ARCTIC
figure(cnt); cnt = cnt+1; 
set(gcf,'Position',[100 100 2000 500]); hold on; 
hfree= plot((SIA_AR_avg_FREE-mean_ARF)./1e12, 'Color', [0.99 0.6 0.1] , 'LineWidth', 5);

x2 = [x, fliplr(x)]; inBetween = [(SIA_AR_min-mean_AR)./1e12, fliplr((SIA_AR_max-mean_AR)./1e12)]; hhh=fill(x2, inBetween, [0.8 0.8 0.8]); set(hhh,'facealpha',.6);
set(hhh,'EdgeColor','none'); 
hreana= plot((SIA_AR_avg-mean_AR)./1e12, 'k' , 'LineWidth', 8);
hobs=plot((SIA_AR_obs_ann-mean(SIA_AR_obs_ann(1,11:end)))./1e12, 'g' , 'LineWidth', 5, 'LineStyle', ':');
%
for iii=1:5
temp1=(SIA_AR_min_hind_anom(records(1,iii)-1,:)./1e12);
temp2=(SIA_AR_max_hind_anom(records(1,iii)-1,:)./1e12);
x3 = [targx(iii,:), fliplr(targx(iii,:))]; inBetween = [temp1, fliplr(temp2)]; hh=fill(x3, inBetween, 'b'); 
set(hh,'facealpha',.2); set(hh,'EdgeColor','none')
end; 
iii=6
temp1=(SIA_AR_min_hind_anom(records(1,iii)-1,1:8)./1e12); temp2=(SIA_AR_max_hind_anom(records(1,iii)-1,1:8)./1e12);
x3 = [targx(iii,1:8), fliplr(targx(iii,1:8))]; inBetween = [temp1, fliplr(temp2)]; hh=fill(x3, inBetween, 'b');
set(hh,'facealpha',.2); set(hh,'EdgeColor','none');
hhind=plot(targx(1,:), squeeze(SIA_AR_avg_hind_anom(records(1,1)-1,:)./1e12), 'b' , 'LineWidth', 8);
for iii=2:5
plot(targx(iii,:), squeeze(SIA_AR_avg_hind_anom(records(1,iii)-1,:)./1e12), 'b' , 'LineWidth', 8);
end;
iii=6; plot(targx(iii,1:8), squeeze(SIA_AR_avg_hind_anom(records(1,iii)-1,1:8)./1e12), 'b' , 'LineWidth', 8);
grid on;
set(gca,'XTick', [0 10 20 30 40 50 60]+1); set(gca,'XTickLabel', [1950 1960 1970 1980 1990 2000 2010]);
set(gca,'fontsize',20); ttitle= ['Arctic SIA ' CASE  ' (drift adj.)']; title(ttitle, 'FontSize',24);
xlabel('years'); ylabel('SIA [mio km^2]'); xlim([11 69]); hleglines = [hfree(1) hreana(1) hhind(1) hobs(1)];
hleg = legend(hleglines,'historical', ['assim-' CASE], ['hindcast-' CASE],'HadISST1');
hleg.Location='southwest'; outputname=[OUTPUTPATH 'SIA_Arctic_driftadj_time' CASE 'F1.eps'];
print(gcf,'-depsc2',outputname);


%SOUTHERN OCEAN
figure(cnt); cnt = cnt+1; set(gcf,'Position',[100 100 2000 500]); hold on; 
hfree= plot((SIA_SO_avg_FREE-mean_SOF)./1e12, 'Color', [0.99 0.6 0.1] , 'LineWidth', 5);

x2 = [x, fliplr(x)]; inBetween = [(SIA_SO_min-mean_SO)./1e12, fliplr((SIA_SO_max-mean_SO)./1e12)]; hhh=fill(x2, inBetween, [0.8 0.8 0.8]); set(hhh,'facealpha',.6); set(hhh,'EdgeColor','none'); 
hreana= plot((SIA_SO_avg-mean_SO)./1e12, 'k' , 'LineWidth', 8);
hobs=plot((SIA_SO_obs_ann-mean(SIA_SO_obs_ann(1,11:end)))./1e12, 'g' , 'LineWidth', 5, 'LineStyle', ':');
%

for iii=1:5
temp1=(SIA_SO_min_hind_anom(records(1,iii)-1,:)./1e12);
temp2=(SIA_SO_max_hind_anom(records(1,iii)-1,:)./1e12);
x3 = [targx(iii,:), fliplr(targx(iii,:))]; inBetween = [temp1, fliplr(temp2)]; hh=fill(x3, inBetween, 'b'); %[0.9+0.1 0.3+0.1 0.9+0.1]);
set(hh,'facealpha',.12); set(hh,'EdgeColor','none')
end; 
iii=6
temp1=(SIA_SO_min_hind_anom(records(1,iii)-1,1:8)./1e12); temp2=(SIA_SO_max_hind_anom(records(1,iii)-1,1:8)./1e12);
x3 = [targx(iii,1:8), fliplr(targx(iii,1:8))]; inBetween = [temp1, fliplr(temp2)]; hh=fill(x3, inBetween, 'b'); 
set(hh,'facealpha',.2); set(hh,'EdgeColor','none');

hhind=plot(targx(1,:), squeeze(SIA_SO_avg_hind_anom(records(1,1)-1,:)./1e12), 'b' , 'LineWidth', 8);
for iii=2:5
plot(targx(iii,:), squeeze(SIA_SO_avg_hind_anom(records(1,iii)-1,:)./1e12), 'b' , 'LineWidth', 8);
end;
iii=6; plot(targx(iii,1:8), squeeze(SIA_SO_avg_hind_anom(records(1,iii)-1,1:8)./1e12), 'b' , 'LineWidth', 8);

grid on;
set(gca,'XTick', [0 10 20 30 40 50 60]+1); set(gca,'XTickLabel', [1950 1960 1970 1980 1990 2000 2010]);
set(gca,'fontsize',20); ttitle= ['Southern Ocean SIA ' CASE  ' (drift adj.)']; title(ttitle, 'FontSize',24);
xlabel('years'); ylabel('SIA [mio km^2]'); xlim([11 69]); hleglines = [hfree(1) hreana(1) hhind(1) hobs(1)];
hleg = legend(hleglines,'historical', ['assim-' CASE], ['hindcast-' CASE],'HadISST1');
hleg.Location='southwest'; outputname=[OUTPUTPATH 'SIA_SOcean_driftadj_time' CASE 'F1.eps'];
print(gcf,'-depsc2',outputname);
%===



end;






