% 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 ='i1';
%reanalysis: 'SIA_Arctic', 'SIA_SOcean', dimensions: (30,nb_yr*12);
load([OUTPUTPATH 'SIV_' CASE '_time.mat']);
%hindcast: 'SIA_Arctic_hind', 'SIA_SOcean_hind', dimensions: (10,nb_yr, 10*12)
load([OUTPUTPATH 'SIV_' CASE 'hind_time.mat']);
% observations 'SIA_HadISST_Arctic', 'SIA_HadISST_SOcean', dimensions: (1,nb_yr*12);
load([OUTPUTPATH 'SIV_PIOMAS_time.mat'])
%free: 'SIV_Arctic_FREE', 'SIV_SOcean_FREE', dimensions: (30,(2014-1950+1)*12);
load([OUTPUTPATH 'SIV_FREE_time.mat']);

%=== annual avgs
%observations and reanalysis start from 1950 jan
%hindcast start from 1960, nov
nmths=828;
nyrs=69;
SIV_SO_ann=zeros(30,69);SIV_AR_ann=zeros(30,69);
SIV_SO_ann_FREE=zeros(30,69);SIV_AR_ann_FREE=zeros(30,69);

SIV_AR_obs_ann=zeros(1,2018-1979+1);
SIV_SO_ann_hind=zeros(10,59,9);SIV_AR_ann_hind=zeros(10,59,9);%start 10 years later in 1960

for i=1:nyrs
  for mem=1:30
    SIV_SO_ann(mem,i)=mean(SIV_SOcean(mem,(i-1)*12+1:i*12));
    SIV_AR_ann(mem,i)=mean(SIV_Arctic(mem,(i-1)*12+1:i*12));
  end; 
end;

for i=1:nyrs
  for mem=1:30
    SIV_SO_ann_FREE(mem,i)=mean(SIV_SOcean_FREE(mem,(i-1)*12+1:i*12));
    SIV_AR_ann_FREE(mem,i)=mean(SIV_Arctic_FREE(mem,(i-1)*12+1:i*12));
  end; 
end;


for i=1:2018-1979+1
  SIV_AR_obs_ann(1,i)=SIV_PIOMAS_Arctic(i,1); 
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;
    SIV_SO_ann_hind(mem,i,lyr)=mean(SIV_SOcean_hind(mem,i,start:ende));
    SIV_AR_ann_hind(mem,i,lyr)=mean(SIV_Arctic_hind(mem,i,start:ende));
   end;
  end; 
end;

%=== determine min, max, avg
SIV_SO_min=zeros(1,69);SIV_AR_min=zeros(1,69);
SIV_SO_max=zeros(1,69);SIV_AR_max=zeros(1,69);
SIV_SO_avg=zeros(1,69);SIV_AR_avg=zeros(1,69);

SIV_SO_min_FREE=zeros(1,69);SIV_AR_min_FREE=zeros(1,69);
SIV_SO_max_FREE=zeros(1,69);SIV_AR_max_FREE=zeros(1,69);
SIV_SO_avg_FREE=zeros(1,69);SIV_AR_avg_FREE=zeros(1,69);


SIV_SO_min_hind=zeros(59,9);SIV_AR_min_hind=zeros(59,9);
SIV_SO_max_hind=zeros(59,9);SIV_AR_max_hind=zeros(59,9);
SIV_SO_avg_hind=zeros(59,9);SIV_AR_avg_hind=zeros(59,9);

for i=1:nyrs
    SIV_SO_min(1,i)=min(SIV_SO_ann(:,i));
    SIV_AR_min(1,i)=min(SIV_AR_ann(:,i));
    SIV_SO_max(1,i)=max(SIV_SO_ann(:,i));
    SIV_AR_max(1,i)=max(SIV_AR_ann(:,i));
    SIV_SO_avg(1,i)=mean(SIV_SO_ann(:,i));
    SIV_AR_avg(1,i)=mean(SIV_AR_ann(:,i));
end;



for i=1:nyrs
    SIV_SO_min_FREE(1,i)=min(SIV_SO_ann_FREE(:,i));
    SIV_AR_min_FREE(1,i)=min(SIV_AR_ann_FREE(:,i));
    SIV_SO_max_FREE(1,i)=max(SIV_SO_ann_FREE(:,i));
    SIV_AR_max_FREE(1,i)=max(SIV_AR_ann_FREE(:,i));
    SIV_SO_avg_FREE(1,i)=mean(SIV_SO_ann_FREE(:,i));
    SIV_AR_avg_FREE(1,i)=mean(SIV_AR_ann_FREE(:,i));
end;


for i=1:nyrs-10
  for ly=1:9
    SIV_SO_min_hind(i,ly)=min(SIV_SO_ann_hind(:,i,ly));
    SIV_AR_min_hind(i,ly)=min(SIV_AR_ann_hind(:,i,ly));
    SIV_SO_max_hind(i,ly)=max(SIV_SO_ann_hind(:,i,ly));
    SIV_AR_max_hind(i,ly)=max(SIV_AR_ann_hind(:,i,ly));
    SIV_SO_avg_hind(i,ly)=mean(SIV_SO_ann_hind(:,i,ly));
    SIV_AR_avg_hind(i,ly)=mean(SIV_AR_ann_hind(:,i,ly));
  end;
end;


%=======
% pick every 10 years and plot min,max,avg hull
porecords=1979-1950+1 :2018-1950+1;
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; 
%x0 = [x, fliplr(x)]; inBetween = [SIV_AR_min_FREE./1e12, fliplr(SIV_AR_max_FREE./1e12)]; h4=fill(x0, inBetween, [0.9 0.6 0.1]);
hfree= plot(SIV_AR_avg_FREE./1e12, 'Color', [0.99 0.6 0.1] , 'LineWidth', 5);
%set(h4,'EdgeColor','none') set(h4,'facealpha',.2); 
x2 = [x, fliplr(x)]; inBetween = [SIV_AR_min./1e12, fliplr(SIV_AR_max./1e12)]; hhh=fill(x2, inBetween, [0.8 0.8 0.8]);
set(hhh,'EdgeColor','none');set(hhh,'facealpha',.6); hreana= plot(SIV_AR_avg./1e12, 'k' , 'LineWidth', 8);
hobs=plot(porecords, SIV_AR_obs_ann, 'Color',[0.0 0.75 0.1] , 'LineWidth', 5, 'LineStyle', ':');
%
for iii=1:5
temp1=(SIV_AR_min_hind(records(1,iii)-1,:)./1e12);
temp2=(SIV_AR_max_hind(records(1,iii)-1,:)./1e12);
x3 = [targx(iii,:), fliplr(targx(iii,:))]; inBetween = [temp1, fliplr(temp2)]; hh=fill(x3, inBetween, 'r'); 
set(hh,'facealpha',.2); set(hh,'EdgeColor','none')
end; 
iii=6
temp1=(SIV_AR_min_hind(records(1,iii)-1,1:8)./1e12); temp2=(SIV_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, 'r');
set(hh,'facealpha',.12); set(hh,'EdgeColor','none');

hhind=plot(targx(1,:), squeeze(SIV_AR_avg_hind(records(1,1)-1,:)./1e12), 'Color', 'r', 'LineWidth', 8);
for iii=2:5
plot(targx(iii,:), squeeze(SIV_AR_avg_hind(records(1,iii)-1,:)./1e12), 'r' , 'LineWidth', 8);
end;
iii=6; plot(targx(iii,1:8), squeeze(SIV_AR_avg_hind(records(1,iii)-1,1:8)./1e12), 'r' , '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 SIV ' CASE ]; title(ttitle, 'FontSize',24);
xlabel('years'); ylabel('SIV [1000 km^3]');
xlim([1 69]); ylim([10 80]);
hleglines = [hfree(1) hreana(1) hhind(1) hobs(1)];
hleg = legend(hleglines,'historical', ['assim-' CASE], ['hindcast-' CASE],'PIOMAS');%,'Location','eastoutside', 'FontSize',18);
hleg.Location='southwest';
outputname=[OUTPUTPATH 'SIV_Arctic_time_PIOMAS' 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(SIV_SO_avg_FREE./1e12, 'Color', [0.99 0.6 0.1] , 'LineWidth', 5);
x2 = [x, fliplr(x)]; inBetween = [SIV_SO_min./1e12, fliplr(SIV_SO_max./1e12)]; hhh=fill(x2, inBetween, [0.8 0.8 0.8]);
set(hhh,'EdgeColor','none'); set(hhh,'facealpha',.6); 
hreana= plot(SIV_SO_avg./1e12, 'k' , 'LineWidth', 8);
%
for iii=1:5
temp1=(SIV_SO_min_hind(records(1,iii)-1,:)./1e12); temp2=(SIV_SO_max_hind(records(1,iii)-1,:)./1e12);
x3 = [targx(iii,:), fliplr(targx(iii,:))]; inBetween = [temp1, fliplr(temp2)]; hh=fill(x3, inBetween, 'r');
set(hh,'facealpha',.12); set(hh,'EdgeColor','none');
end; 
iii=6
temp1=(SIV_SO_min_hind(records(1,iii)-1,1:8)./1e12); temp2=(SIV_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, 'r'); 
set(hh,'facealpha',.2); set(hh,'EdgeColor','none');


hhind=plot(targx(1,:), squeeze(SIV_SO_avg_hind(records(1,1)-1,:)./1e12), 'r' , 'LineWidth', 8);
for iii=2:5
plot(targx(iii,:), squeeze(SIV_SO_avg_hind(records(1,iii)-1,:)./1e12), 'r' , 'LineWidth', 8);
end;
iii=6; plot(targx(iii,1:8), squeeze(SIV_SO_avg_hind(records(1,iii)-1,1:8)./1e12), 'r' , '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 SIV ' CASE ]; title(ttitle, 'FontSize',24);
xlabel('years'); ylabel('SIV [1000 km^3]');
xlim([1 69]); ylim([10 30])
hleglines = [hfree(1) hreana(1) hhind(1) ];
hleg = legend(hleglines,'historical', ['assim-' CASE], ['hindcast-' CASE]);%,'Location','eastoutside', 'FontSize',18);
hleg.Location='southwest';
outputname=[OUTPUTPATH 'SIV_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)
%
% respective time period where climatlogy is taken from: 1979-2018 (PIOMAS)
%   for hindcast, 1979 is entry number 20
%   for reanalyis, 1979 is entry number 30
clim_t_AR=zeros(1,9);
clim_t_SO=zeros(1,9);
for ly=1:9
  clim_t_AR(1,ly) = mean(SIV_AR_avg_hind(20:end,ly));  % for PIOMAS!
  clim_t_SO(1,ly) = mean(SIV_SO_avg_hind(20:end,ly));  % for PIOMAS!
end;
SIV_SO_avg_hind_anom=zeros(59,9); 
SIV_AR_avg_hind_anom=zeros(59,9);
SIV_SO_min_hind_anom=zeros(59,9); 
SIV_AR_min_hind_anom=zeros(59,9);
SIV_SO_max_hind_anom=zeros(59,9); 
SIV_AR_max_hind_anom=zeros(59,9);
for i=1:nyrs-10
for ly=1:9
  SIV_SO_avg_hind_anom(i,ly)=SIV_SO_avg_hind(i,ly)-clim_t_SO(1,ly);
  SIV_AR_avg_hind_anom(i,ly)=SIV_AR_avg_hind(i,ly)-clim_t_AR(1,ly);
  SIV_SO_min_hind_anom(i,ly)=SIV_SO_min_hind(i,ly)-clim_t_SO(1,ly);
  SIV_AR_min_hind_anom(i,ly)=SIV_AR_min_hind(i,ly)-clim_t_AR(1,ly);
  SIV_SO_max_hind_anom(i,ly)=SIV_SO_max_hind(i,ly)-clim_t_SO(1,ly);
  SIV_AR_max_hind_anom(i,ly)=SIV_AR_max_hind(i,ly)-clim_t_AR(1,ly);

end; end; 
mean_AR=mean(SIV_AR_avg(1,30:end)); %over same period as hindcast
mean_SO=mean(SIV_SO_avg(1,30:end)); %over same period as hindcast
mean_ARF=mean(SIV_AR_avg_FREE(1,30:end)); %over same period as hindcast
mean_SOF=mean(SIV_SO_avg_FREE(1,30:end)); %over same period as hindcast

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

x2 = [x, fliplr(x)]; inBetween = [(SIV_AR_min-mean_AR)./1e12, fliplr((SIV_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((SIV_AR_avg-mean_AR)./1e12, 'k' , 'LineWidth', 8);
hobs=plot(porecords,(SIV_AR_obs_ann-mean(SIV_AR_obs_ann(1,:))), 'Color',[0.0 0.75 0.1] , 'LineWidth', 5, 'LIneStyle', ':');
%
for iii=1:5
temp1=(SIV_AR_min_hind_anom(records(1,iii)-1,:)./1e12); temp2=(SIV_AR_max_hind_anom(records(1,iii)-1,:)./1e12);
x3 = [targx(iii,:), fliplr(targx(iii,:))]; inBetween = [temp1, fliplr(temp2)]; hh=fill(x3, inBetween, 'r'); %[0.9+0.1 0.3+0.1 0.9+0.1]);
set(hh,'facealpha',.12); set(hh,'EdgeColor','none')
end; 
iii=6
temp1=(SIV_AR_min_hind_anom(records(1,iii)-1,1:8)./1e12); temp2=(SIV_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, 'r'); %[0.9+0.1 0.3+0.1 0.9+0.1]);
set(hh,'facealpha',.2); set(hh,'EdgeColor','none');

hhind=plot(targx(1,:), squeeze(SIV_AR_avg_hind_anom(records(1,1)-1,:)./1e12), 'r' , 'LineWidth', 8);
for iii=2:5
plot(targx(iii,:), squeeze(SIV_AR_avg_hind_anom(records(1,iii)-1,:)./1e12), 'r' , 'LineWidth', 8);
end;
iii=6; plot(targx(iii,1:8), squeeze(SIV_AR_avg_hind_anom(records(1,iii)-1,1:8)./1e12), 'r' , '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 SIV ' CASE  ' (drift adj.)']; title(ttitle, 'FontSize',24);
xlabel('years'); ylabel('SIV [1000 km^3]'); xlim([11 69]); ylim([-15 20]); 
hleglines = [hfree(1) hreana(1) hhind(1) hobs(1)];
hleg = legend(hleglines,'historical', ['assim-' CASE], ['hindcast-' CASE],'PIOMAS');%,'Location','eastoutside', 'FontSize',18);
hleg.Location='southwest';
outputname=[OUTPUTPATH 'SIV_Arctic_driftadj_PIOMAS_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((SIV_SO_avg_FREE-mean_SOF)./1e12, 'Color', [0.99 0.6 0.1] , 'LineWidth', 5);
x2 = [x, fliplr(x)]; inBetween = [(SIV_SO_min-mean_SO)./1e12, fliplr((SIV_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((SIV_SO_avg-mean_SO)./1e12, 'k' , 'LineWidth', 8);
%
for iii=1:5
temp1=(SIV_SO_min_hind_anom(records(1,iii)-1,:)./1e12);  temp2=(SIV_SO_max_hind_anom(records(1,iii)-1,:)./1e12);
x3 = [targx(iii,:), fliplr(targx(iii,:))]; inBetween = [temp1, fliplr(temp2)]; hh=fill(x3, inBetween, 'r'); 
set(hh,'facealpha',.12); set(hh,'EdgeColor','none')
end; 
iii=6
temp1=(SIV_SO_min_hind_anom(records(1,iii)-1,1:8)./1e12); temp2=(SIV_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, 'r'); 
set(hh,'facealpha',.2); set(hh,'EdgeColor','none');

hhind=plot(targx(1,:), squeeze(SIV_SO_avg_hind_anom(records(1,1)-1,:)./1e12), 'r' , 'LineWidth', 8);
for iii=2:5
plot(targx(iii,:), squeeze(SIV_SO_avg_hind_anom(records(1,iii)-1,:)./1e12), 'r' , 'LineWidth', 8);
end;
iii=6; plot(targx(iii,1:8), squeeze(SIV_SO_avg_hind_anom(records(1,iii)-1,1:8)./1e12), 'r' , '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 SIV ' CASE  ' (drift adj.)']; title(ttitle, 'FontSize',24);
xlabel('years'); ylabel('SIV [1000 km^3]'); xlim([11 69]);  ylim([-4 10]); 
hleglines = [hfree(1) hreana(1) hhind(1) ];
hleg = legend(hleglines,'historical', ['assim-' CASE], ['hindcast-' CASE]);%,'Location','eastoutside', 'FontSize',18);
hleg.Location='northeast'; outputname=[OUTPUTPATH 'SIV_SOcean_driftadj_time' CASE 'F1.eps'];
print(gcf,'-depsc2',outputname);
%===



end;