%Comparison of energy balance fluxes etc. between standard run and Ta +2

%Note first run validation code on climate change outputs then run through
%the T_and_C_Peru_Analysis code and save everything


%Load data
%Both updated 23/02/2021 to account for change to vapour pressure and dew
%point
SD=load('SD_all_2.mat'); %Note this is the climate change standard run which is not the same as the all other analysis final run due to differences in albedo and ice depth
Ta2=load('Ta2_all_2.mat'); %Cool that loads everything into a structure %Ta plus 2C no Pr change

%Staked bar chart of seasonal fluxes

%First sort out data
All_A_fluxes_Ta2_diff=Ta2.All_A_fluxes-SD.All_A_fluxes; %Scenario - standard run
%Restructure so goes standard, scenario, difference then next station
Gap=zeros(1,8);
All_SD_Ta2_A_fluxes=[SD.All_A_fluxes(1:4,:);Ta2.All_A_fluxes(1:4,:);All_A_fluxes_Ta2_diff(1:4,:);SD.All_A_fluxes(5:8,:);Ta2.All_A_fluxes(5:8,:);All_A_fluxes_Ta2_diff(5:8,:);SD.All_A_fluxes(9:12,:);Ta2.All_A_fluxes(9:12,:);All_A_fluxes_Ta2_diff(9:12,:);SD.All_A_fluxes(13:16,:);Ta2.All_A_fluxes(13:16,:);All_A_fluxes_Ta2_diff(13:16,:);SD.All_A_fluxes(17:19,:);Gap;Ta2.All_A_fluxes(17:19,:);Gap;All_A_fluxes_Ta2_diff(17:19,:)];
All_SD_Ta2_A_fluxes(52,8)=NaN; All_SD_Ta2_A_fluxes(56,8)=NaN; 

%Right so want seasonal ablation types stacked and the differences.
%First make a table from A_All

%So this is the standard run one
SD.SG_A_all_T=array2table(SD.SG_A_all,'VariableNames',SD.SG_T_labels); SD.SG_A_all_T.Sub=SD.SG_A_all_T.EICE+SD.SG_A_all_T.ESN; %Calculate total sublimation mean
SD.AG_A_all_T=array2table(SD.AG_A_all,'VariableNames',SD.AG_T_labels); SD.AG_A_all_T.Sub=SD.AG_A_all_T.EICE+SD.AG_A_all_T.ESN; %Calculate total sublimation mean
SD.CG_A_all_T=array2table(SD.CG_A_all,'VariableNames',SD.CG_T_labels); SD.CG_A_all_T.Sub=SD.CG_A_all_T.EICE+SD.CG_A_all_T.ESN; %Calculate total sublimation mean
SD.QQG_A_all_T=array2table(SD.QQG_A_all,'VariableNames',SD.QQG_T_labels); SD.QQG_A_all_T.Sub=SD.QQG_A_all_T.EICE+SD.QQG_A_all_T.ESN; %Calculate total sublimation mean
SD.QIC_A_all_T=array2table(SD.QIC_A_all,'VariableNames',SD.QIC_T_labels); SD.QIC_A_all_T.Sub=SD.QIC_A_all_T.EICE+SD.QIC_A_all_T.ESN; %Calculate total sublimation mean

SD.SG_A_ablat_type(1,:)=SD.SG_A_all_T(:,{'Imelt','Smelt','Sub'});
SD.SG_A_ablat_type(2,:)=SD.SG_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Dry season
SD.SG_A_ablat_type(3,:)=SD.SG_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Wet season
SD.AG_A_ablat_type(1,:)=SD.AG_A_all_T(:,{'Imelt','Smelt','Sub'});
SD.AG_A_ablat_type(2,:)=SD.AG_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Dry season
SD.AG_A_ablat_type(3,:)=SD.AG_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Wet season
SD.CG_A_ablat_type(1,:)=SD.CG_A_all_T(:,{'Imelt','Smelt','Sub'});
SD.CG_A_ablat_type(2,:)=SD.CG_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Dry season
SD.CG_A_ablat_type(3,:)=SD.CG_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Wet season
SD.QIC_A_ablat_type(1,:)=SD.QIC_A_all_T(:,{'Imelt','Smelt','Sub'});
SD.QIC_A_ablat_type(2,:)=SD.QIC_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Dry season
SD.QIC_A_ablat_type(3,:)=SD.QIC_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Wet season
SD.QQG_A_ablat_type(1,:)=SD.QQG_A_all_T(:,{'Imelt','Smelt','Sub'});
SD.QQG_A_ablat_type(2,:)=SD.QQG_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Dry season
SD.QQG_A_ablat_type(3,:)=SD.QQG_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Wet season

Gap=zeros(1,3);
SD.All_A_ablat_type=[table2array(SD.SG_A_ablat_type);Gap;table2array(SD.AG_A_ablat_type);Gap;table2array(SD.CG_A_ablat_type);Gap;table2array(SD.QQG_A_ablat_type);Gap;table2array(SD.QIC_A_ablat_type)];

%So this is for Ta + 2
Ta2.SG_A_all_T=array2table(Ta2.SG_A_all,'VariableNames',Ta2.SG_T_labels); Ta2.SG_A_all_T.Sub=Ta2.SG_A_all_T.EICE+Ta2.SG_A_all_T.ESN; %Calculate total sublimation mean
Ta2.AG_A_all_T=array2table(Ta2.AG_A_all,'VariableNames',Ta2.AG_T_labels); Ta2.AG_A_all_T.Sub=Ta2.AG_A_all_T.EICE+Ta2.AG_A_all_T.ESN; %Calculate total sublimation mean
Ta2.CG_A_all_T=array2table(Ta2.CG_A_all,'VariableNames',Ta2.CG_T_labels); Ta2.CG_A_all_T.Sub=Ta2.CG_A_all_T.EICE+Ta2.CG_A_all_T.ESN; %Calculate total sublimation mean
Ta2.QQG_A_all_T=array2table(Ta2.QQG_A_all,'VariableNames',Ta2.QQG_T_labels); Ta2.QQG_A_all_T.Sub=Ta2.QQG_A_all_T.EICE+Ta2.QQG_A_all_T.ESN; %Calculate total sublimation mean
Ta2.QIC_A_all_T=array2table(Ta2.QIC_A_all,'VariableNames',Ta2.QIC_T_labels); Ta2.QIC_A_all_T.Sub=Ta2.QIC_A_all_T.EICE+Ta2.QIC_A_all_T.ESN; %Calculate total sublimation mean

Ta2.SG_A_ablat_type(1,:)=Ta2.SG_A_all_T(:,{'Imelt','Smelt','Sub'});
Ta2.SG_A_ablat_type(2,:)=Ta2.SG_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Dry season
Ta2.SG_A_ablat_type(3,:)=Ta2.SG_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Wet season
Ta2.AG_A_ablat_type(1,:)=Ta2.AG_A_all_T(:,{'Imelt','Smelt','Sub'});
Ta2.AG_A_ablat_type(2,:)=Ta2.AG_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Dry season
Ta2.AG_A_ablat_type(3,:)=Ta2.AG_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Wet season
Ta2.CG_A_ablat_type(1,:)=Ta2.CG_A_all_T(:,{'Imelt','Smelt','Sub'});
Ta2.CG_A_ablat_type(2,:)=Ta2.CG_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Dry season
Ta2.CG_A_ablat_type(3,:)=Ta2.CG_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Wet season
Ta2.QIC_A_ablat_type(1,:)=Ta2.QIC_A_all_T(:,{'Imelt','Smelt','Sub'});
Ta2.QIC_A_ablat_type(2,:)=Ta2.QIC_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Dry season
Ta2.QIC_A_ablat_type(3,:)=Ta2.QIC_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Wet season
Ta2.QQG_A_ablat_type(1,:)=Ta2.QQG_A_all_T(:,{'Imelt','Smelt','Sub'});
Ta2.QQG_A_ablat_type(2,:)=Ta2.QQG_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Dry season
Ta2.QQG_A_ablat_type(3,:)=Ta2.QQG_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub'}); %Wet season

Gap=zeros(1,3);
Ta2.All_A_ablat_type=[table2array(Ta2.SG_A_ablat_type);Gap;table2array(Ta2.AG_A_ablat_type);Gap;table2array(Ta2.CG_A_ablat_type);Gap;table2array(Ta2.QQG_A_ablat_type);Gap;table2array(Ta2.QIC_A_ablat_type)];

%Then find the difference
All_A_ablat_type_Ta2_diff=Ta2.All_A_ablat_type-SD.All_A_ablat_type; %Scenario - standard run

%Figures
figure5=figure();
%First difference in fluxes
ax=subplot(1,2,1);
All_labels=categorical({'A','D','W','A','D','W','A','D','W','A','D','W','A','D','W'});
ax.ColorOrder = SD.Fluxes_colours; hold on; bar(All_A_fluxes_Ta2_diff(:,1:7),'stacked'); hold on; %Not including dQ
scatter([1:19],All_A_fluxes_Ta2_diff(:,8),22,'w','d','filled','MarkerEdgeColor',SD.Fluxes_colours(1,:));
legend('S^\ast','L^\ast','H','\lambdaE','Q_v','Q_f_m','G','dQ','Location','SouthOutside','Orientation','horizontal','FontSize',11,'NumColumns',4);
xticks([1 2 3 5 6 7 9 10 11 13 14 15 17 18 19]); %Specifing the x ticks with gaps where we don't want a label
xticklabels(All_labels);
ylabel('Difference in energy fluxes (W m^-^2)','FontSize',14); 
ylim ([-10 70]);

%The difference in ablation
ax=subplot(1,2,2);
ax.ColorOrder = SD.Ablat_colours; hold on; bar(All_A_ablat_type_Ta2_diff,'stacked');hold on;
legend('Ice melt','Snow melt','Sublimation','Location','SouthOutside','Orientation','Horizontal','FontSize',11);
xticks([1 2 3 5 6 7 9 10 11 13 14 15 17 18 19]); %Specifing the x ticks with gaps where we don't want a label
xticklabels(All_labels);
ylabel('Difference in mean ablation (mm w.e. h^-^1)','FontSize',14);
ylim([-0.4 1]);

% Create textbox
annotation(figure5,'textbox',[0.13308161190708 0.865740226908373 0.0638967551622418 0.0534591194968551],'String','SG','LineStyle','none','HorizontalAlignment','center','FontSize',11,'FitBoxToText','off');
annotation(figure5,'textbox',[0.199099871410201 0.868763102725367 0.0638967551622417 0.0534591194968551],'String','AG','LineStyle','none','HorizontalAlignment','center','FontSize',11,'FitBoxToText','off');
annotation(figure5,'textbox',[0.265834820284776 0.867643975829325 0.0638967551622417 0.0534591194968551],'String','CG','LineStyle','none','HorizontalAlignment','center','FontSize',11,'FitBoxToText','off');
annotation(figure5,'textbox',[0.331022775904354 0.870625231224564 0.0638967551622421 0.0534591194968551],'String','QQG','LineStyle','none','HorizontalAlignment','center','FontSize',11,'FitBoxToText','off');
annotation(figure5,'textbox',[0.395353466200315 0.868686027870267 0.0638967551622417 0.0534591194968551],'String','QIC','LineStyle','none','HorizontalAlignment','center','FontSize',11,'FitBoxToText','off');
annotation(figure5,'textbox',[0.573064209669283 0.864433037365889 0.0638967551622419 0.0534591194968551],'String','SG','LineStyle','none','HorizontalAlignment','center','FontSize',11,'FitBoxToText','off');
annotation(figure5,'textbox',[0.638822912524148 0.869786656801085 0.0638967551622414 0.0534591194968551],'String','AG','LineStyle','none','HorizontalAlignment','center','FontSize',11,'FitBoxToText','off');
annotation(figure5,'textbox',[0.705694501263937 0.867643975829325 0.0638967551622417 0.0534591194968551],'String','CG','LineStyle','none','HorizontalAlignment','center','FontSize',11,'FitBoxToText','off');
annotation(figure5,'textbox',[0.770060982584533 0.871932420767048 0.0638967551622419 0.0534591194968551],'String','QQG','LineStyle','none','HorizontalAlignment','center','FontSize',11,'FitBoxToText','off');
annotation(figure5,'textbox',[0.834856096566939 0.868686027870267 0.0638967551622416 0.0534591194968551],'String','QIC','LineStyle','none','HorizontalAlignment','center','FontSize',11,'FitBoxToText','off');

annotation(figure5,'textbox',[0.0772966138019717 0.903267973856209 0.019145735105015 0.0705882352941175],'String',{'a'},'FontSize',16,'FitBoxToText','off','EdgeColor','none');
annotation(figure5,'textbox',[0.516645092156022 0.90718954248366 0.019145735105015 0.0705882352941175],'String','b','FontSize',16,'FitBoxToText','off','EdgeColor','none');

%Looking at differences in precipitation phase
Pr_perc_sno=NaN(5,3);
Pr_perc_sno(1,1)=(SD.SG_S_all_T.Pr_sno/(SD.SG_S_all_T.Pr_liq+SD.SG_S_all_T.Pr_sno))*100; %SG SD
Pr_perc_sno(1,2)=(Ta2.SG_S_all_T.Pr_sno/(Ta2.SG_S_all_T.Pr_liq+Ta2.SG_S_all_T.Pr_sno))*100; %SG Ta2
Pr_perc_sno(1,3)=Pr_perc_sno(1,2)-Pr_perc_sno(1,1); %Difference scenario - standard
Pr_perc_sno(2,1)=(SD.AG_S_all_T.Pr_sno/(SD.AG_S_all_T.Pr_liq+SD.AG_S_all_T.Pr_sno))*100; %AG SD
Pr_perc_sno(2,2)=(Ta2.AG_S_all_T.Pr_sno/(Ta2.AG_S_all_T.Pr_liq+Ta2.AG_S_all_T.Pr_sno))*100; %AG Ta2
Pr_perc_sno(2,3)=Pr_perc_sno(2,2)-Pr_perc_sno(2,1); %Difference scenario - standard
Pr_perc_sno(3,1)=(SD.CG_S_all_T.Pr_sno/(SD.CG_S_all_T.Pr_liq+SD.CG_S_all_T.Pr_sno))*100; %CG SD
Pr_perc_sno(3,2)=(Ta2.CG_S_all_T.Pr_sno/(Ta2.CG_S_all_T.Pr_liq+Ta2.CG_S_all_T.Pr_sno))*100; %CG Ta2
Pr_perc_sno(3,3)=Pr_perc_sno(3,2)-Pr_perc_sno(3,1); %Difference scenario - standard
Pr_perc_sno(4,1)=(SD.QQG_S_all_T.Pr_sno/(SD.QQG_S_all_T.Pr_liq+SD.QQG_S_all_T.Pr_sno))*100; %QQG SD
Pr_perc_sno(4,2)=(Ta2.QQG_S_all_T.Pr_sno/(Ta2.QQG_S_all_T.Pr_liq+Ta2.QQG_S_all_T.Pr_sno))*100; %QQG Ta2
Pr_perc_sno(4,3)=Pr_perc_sno(4,2)-Pr_perc_sno(4,1); %Difference scenario - standard
Pr_perc_sno(5,1)=(SD.QIC_S_all_T.Pr_sno/(SD.QIC_S_all_T.Pr_liq+SD.QIC_S_all_T.Pr_sno))*100; %QIC SD
Pr_perc_sno(5,2)=(Ta2.QIC_S_all_T.Pr_sno/(Ta2.QIC_S_all_T.Pr_liq+Ta2.QIC_S_all_T.Pr_sno))*100; %QIC Ta2
Pr_perc_sno(5,3)=Pr_perc_sno(5,2)-Pr_perc_sno(5,1); %Difference scenario - standard

%looking at albedo difference
Alb_mean_all=NaN(5,3);
Alb_mean_all(1,1)=nanmean(SD.SG_T.ALB);
Alb_mean_all(1,2)=nanmean(Ta2.SG_T.ALB);
Alb_mean_all(1,3)=Alb_mean_all(1,2)-Alb_mean_all(1,1); %Difference scenario - standard
Alb_mean_all(2,1)=nanmean(SD.AG_T.ALB);
Alb_mean_all(2,2)=nanmean(Ta2.AG_T.ALB);
Alb_mean_all(2,3)=Alb_mean_all(2,2)-Alb_mean_all(2,1); %Difference scenario - standard
Alb_mean_all(3,1)=nanmean(SD.CG_T.ALB);
Alb_mean_all(3,2)=nanmean(Ta2.CG_T.ALB);
Alb_mean_all(3,3)=Alb_mean_all(3,2)-Alb_mean_all(3,1); %Difference scenario - standard
Alb_mean_all(4,1)=nanmean(SD.QQG_T.ALB);
Alb_mean_all(4,2)=nanmean(Ta2.QQG_T.ALB);
Alb_mean_all(4,3)=Alb_mean_all(4,2)-Alb_mean_all(4,1); %Difference scenario - standard
Alb_mean_all(5,1)=nanmean(SD.QIC_T.ALB);
Alb_mean_all(5,2)=nanmean(Ta2.QIC_T.ALB);
Alb_mean_all(5,3)=Alb_mean_all(5,2)-Alb_mean_all(5,1); %Difference scenario - standard


