% initial development : Madlen Kimmritz (madlen.kimmritz@nersc.no), Nov 2019
% Plot time development of SIA (annual avgs) for
%  - reanalysis (avg, max, min)
%  - free
%  - observation
%  - hindcast (avg,max,min) branched every n years (n=10)
%  - optional: drift adjustment
% for southern and northern hemisphere
% As observations for SIT are provided only for months Oct-April, 
% we consider annual values that are averages over Oct-Apr
%===
% period of available observations winter01: 2002/03 - winter15: 2016/17
% winter01: 1950/51, winter68: 2017/2018
%====
addpath('/cluster/home/madimm/matlab')
%==time for reanalysis
yr_start=1951; yr_end  =2018; nyrs=yr_end-yr_start + 1;

%=== read data
OUTPUTPATH='/cluster/home/madimm/PROGRAMS_PaperIngoDCPP/PICS/';
CASE ='i1';
load([OUTPUTPATH 'SIV_' CASE '_time.mat']);
load([OUTPUTPATH 'SIV_FREE_time.mat']);
load([OUTPUTPATH 'SIV_' CASE 'hind_time.mat']);
load([OUTPUTPATH 'SIV_OBSmerged_time.mat'])

%=== annual avgs= avg over Oct yr1 - Apr yr1+1
%observations and reanalysis start from 1950 jan
%hindcast start from 1960, nov
nmths=828;
SIV_SO_ann=zeros(30,nyrs);SIV_AR_ann=zeros(30,nyrs);
SIV_SO_ann_FREE=zeros(30,nyrs);SIV_AR_ann_FREE=zeros(30,nyrs);
SIV_AR_obs_ann=zeros(1,15);SIV_AR_obs_ann_had=zeros(1,15);
SIV_SO_ann_hind=zeros(10,nyrs-10,9);SIV_AR_ann_hind=zeros(10,nyrs-10,9);

for i=1:nyrs
  for mem=1:30
    start=(i-1)*12+10; ende=start+6;
    SIV_SO_ann(mem,i)=mean(SIV_SOcean(mem,start:ende));
    SIV_AR_ann(mem,i)=mean(SIV_Arctic(mem,start:ende));
  end; 
end;
for i=1:nyrs
  for mem=1:30
    start=(i-1)*12+10; ende=start+6;
    SIV_SO_ann_FREE(mem,i)=mean(SIV_SOcean_FREE(mem,start:ende));
    SIV_AR_ann_FREE(mem,i)=mean(SIV_Arctic_FREE(mem,start:ende));
  end; 
end;

time=0;
%SIV_obsmerged_Arctic_time_merged; SIV_obsmerged_Arctic_time_HadISST;

for i=2002:2016
  time=time+1;
  start=(time-1)*12+10; 
  ende=start+6;
  SIV_AR_obs_ann(1,time)=mean(SIV_obsmerged_Arctic_time_merged(1,start:ende)); 
  %SIV_AR_obs_ann_had(1,time)=mean(SIV_obsmerged_Arctic_time_HadISST(1,start:ende));
end;


for i=1:nyrs-10 
  %hindcast starts 10 years later !=====
  for mem=1:10
   %first and last lead year: special: ly1: avg over nov-apr
   lyr=1; 
   start=1; ende=start+5;
   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));
   %6:apr, 7M,8J,9J,10A,11S,12O
   for lyr=2:9
    start=(lyr-2)*12+12; 
    ende=start+6;
    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,nyrs);SIV_AR_min=zeros(1,nyrs);
SIV_SO_max=zeros(1,nyrs);SIV_AR_max=zeros(1,nyrs);
SIV_SO_avg=zeros(1,nyrs);SIV_AR_avg=zeros(1,nyrs);

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

SIV_SO_min_hind=zeros(nyrs-10,9);SIV_AR_min_hind=zeros(nyrs-10,9);
SIV_SO_max_hind=zeros(nyrs-10,9);SIV_AR_max_hind=zeros(nyrs-10,9);
SIV_SO_avg_hind=zeros(nyrs-10,9);SIV_AR_avg_hind=zeros(nyrs-10,9);

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
    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-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=2003-1951+1 :2017-1951+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; 
hfree= plot(SIV_AR_avg_FREE./1e12, 'Color', [0.99 0.6 0.1] , 'LineWidth', 5);

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./1e12, 'Color',[0.0 0.75 0.1] , 'LineWidth', 5, 'LineStyle', ':');
%
hobs1=plot(porecords, SIV_AR_obs_ann_had./1e12, 'Color',[0.9 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'); %[0.9+0.1 0.3+0.1 0.9+0.1]);
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'); %[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(records(1,1)-1,:)./1e12), 'Color', 'r', 'LineWidth', 5);
for iii=2:5
plot(targx(iii,:), squeeze(SIV_AR_avg_hind(records(1,iii)-1,:)./1e12), 'r' , 'LineWidth', 5);
end;
iii=6; plot(targx(iii,1:8), squeeze(SIV_AR_avg_hind(records(1,iii)-1,1:8)./1e12), 'r' , 'LineWidth', 5);

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',18);
ttitle= ['Arctic SIV dcppA-assim-' CASE ]; title(ttitle, 'FontSize',24);
xlabel('years'); ylabel('SIV [1000 km^3]');
xlim([1 69]); ylim([10 80]);
hleglines = [hreana(1) hhind(1) hobs(1) hobs1(1)];
hleg = legend(hleglines,['dcppA-assim-' CASE],'hindcast','Cryo', 'Cryo+HadISST');%,'Location','eastoutside', 'FontSize',18);
hleg.Location='eastoutside';
outputname=[OUTPUTPATH 'SIV_Arctic_time_merged' CASE '.eps'];
print(gcf,'-depsc2',outputname);


%SOUTHERN OCEAN
figure(cnt); cnt = cnt+1; 
set(gcf,'Position',[100 100 2000 500])
hold on; 
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')
hreana= plot(SIV_SO_avg./1e12, 'k' , 'LineWidth', 5);
%

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'); %[0.9+0.1 0.3+0.1 0.9+0.1]);
set(hh,'facealpha',.2); 
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'); %[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_SO_avg_hind(records(1,1)-1,:)./1e12), 'r' , 'LineWidth', 5);
for iii=2:5
plot(targx(iii,:), squeeze(SIV_SO_avg_hind(records(1,iii)-1,:)./1e12), 'r' , 'LineWidth', 5);
end;
iii=6; plot(targx(iii,1:8), squeeze(SIV_SO_avg_hind(records(1,iii)-1,1:8)./1e12), 'r' , 'LineWidth', 5);

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',18);
ttitle= ['Southern Ocean SIV dcppA-assim-' CASE ]; title(ttitle, 'FontSize',24);
xlabel('years'); ylabel('SIV [1000 km^3]');
xlim([1 69]); ylim([10 30])
hleglines = [hreana(1) hhind(1)];
hleg = legend(hleglines,['dcppA-assim-' CASE],'hindcast');%,'Location','eastoutside', 'FontSize',18);
hleg.Location='eastoutside';
outputname=[OUTPUTPATH 'SIV_SouthOcean_time' CASE '.eps'];
print(gcf,'-depsc2',outputname);



%============================= TODO

driftadj=0;

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

%===========
%=== plot
%ARCTIC
figure(cnt); cnt = cnt+1; 
set(gcf,'Position',[100 100 2000 500]); hold on; 
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,'EdgeColor','none'); 
hreana= plot((SIV_AR_avg-mean_AR)./1e12, 'k' , 'LineWidth', 5);
hobs=plot(porecords,(SIV_AR_obs_ann-mean(SIV_AR_obs_ann(1,:))), 'Color',[0.0 0.75 0.1] , 'LineWidth', 4);
%

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',.2); 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', 5);
for iii=2:5
plot(targx(iii,:), squeeze(SIV_AR_avg_hind_anom(records(1,iii)-1,:)./1e12), 'r' , 'LineWidth', 5);
end;
iii=6; plot(targx(iii,1:8), squeeze(SIV_AR_avg_hind_anom(records(1,iii)-1,1:8)./1e12), 'r' , 'LineWidth', 5);

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',18); ttitle= ['Arctic SIV dcppA-assim-' CASE  ' (drift adj.)']; title(ttitle, 'FontSize',24);
xlabel('years'); ylabel('SIV [1000 km^3]'); xlim([11 69]); ylim([-15 20]); hleglines = [hreana(1) hhind(1) hobs(1)];
hleg = legend(hleglines,['dcppA-assim-' CASE],'hindcast','PIOMAS');%,'Location','eastoutside', 'FontSize',18);
hleg.Location='eastoutside'; outputname=[OUTPUTPATH 'SIV_Arctic_driftadj_PIOMAS_time' CASE '.eps'];
print(gcf,'-depsc2',outputname);


%SOUTHERN OCEAN
figure(cnt); cnt = cnt+1; 
set(gcf,'Position',[100 100 2000 500]); hold on; 
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,'EdgeColor','none'); 
hreana= plot((SIV_SO_avg-mean_SO)./1e12, 'k' , 'LineWidth', 5);
%

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'); %[0.9+0.1 0.3+0.1 0.9+0.1]);
set(hh,'facealpha',.2); 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'); %[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_SO_avg_hind_anom(records(1,1)-1,:)./1e12), 'r' , 'LineWidth', 5);
for iii=2:5
plot(targx(iii,:), squeeze(SIV_SO_avg_hind_anom(records(1,iii)-1,:)./1e12), 'r' , 'LineWidth', 5);
end;
iii=6; plot(targx(iii,1:8), squeeze(SIV_SO_avg_hind_anom(records(1,iii)-1,1:8)./1e12), 'r' , 'LineWidth', 5);

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',18); ttitle= ['Southern Ocean SIV dcppA-assim-' CASE  ' (drift adj.)']; title(ttitle, 'FontSize',24);
xlabel('years'); ylabel('SIV [1000 km^3]'); xlim([11 69]);  ylim([-4 10]); hleglines = [hreana(1) hhind(1) hobs(1)];
hleg = legend(hleglines,['dcppA-assim-' CASE],'hindcast');%,'Location','eastoutside', 'FontSize',18);
hleg.Location='eastoutside'; outputname=[OUTPUTPATH 'SIV_SOcean_driftadj_time' CASE '.eps'];
print(gcf,'-depsc2',outputname);
%===



end;






