%Script to analyse the effect of Ta and Pr variations on T and C outputs
%C Fyffe 25/11/2020 %Edited by C Fyffe 23/08/2021 for depositing into EIDC

%Look at mean ice/snow melt and ablation differences in wet and dry season
%and total across all sites

%Load data from climate runs
load('Climate_Sensitivity_vars_2_ea_change.mat');

%===========================================
%Shallap
%===========================================

%Useful things
SG_datasize=size(SG(1).run.Date,1);

%Create matrix for each variable
SG_vars = [SG(1:35).run];
SG_Smelt=[SG_vars(1:35).Smelt]; %yes! this pulls out both Smelts into a table
SG_Imelt=[SG_vars(1:35).Imelt];
SG_Tmelt=SG_Smelt+SG_Imelt; %Good this works - total melt
SG_ESN=[SG_vars(1:35).ESN];
SG_EICE=[SG_vars(1:35).EICE];
SG_Sub=SG_ESN+SG_EICE; %Sublimation
SG_Abl=SG_Smelt+SG_Imelt+SG_ESN+SG_EICE; %Ablation
SG_Ta=[SG_vars(1:35).Ta];
SG_Pr=[SG_vars(1:35).Pr];
SG_Pr_liq=[SG_vars(1:35).Pr_liq];
SG_Pr_sno=[SG_vars(1:35).Pr_sno];
SG_ICE=[SG_vars(1:35).ICE];
SG_SWE=[SG_vars(1:35).SWE];
SG_ISWE=SG_ICE+SG_SWE;

SG_T=[SG_Smelt SG_Imelt SG_Tmelt SG_ESN SG_EICE SG_Sub SG_Abl SG_Ta SG_Pr SG_Pr_liq SG_Pr_sno SG_ICE SG_SWE SG_ISWE]; %Combine data into one matrix

%Create big table
SG_T=array2table(SG_T,'VariableNames',Table_names_long); %Note use of array2table when turning a matrix into a table
SG_DateTime=datetime(SG_vars(1).Date,'ConvertFrom','datenum');
SG_T = addvars(SG_T,SG_DateTime,'Before',1);
SG_T.Year=year(SG_T.SG_DateTime); %Calculate year
SG_T.Month=month(SG_T.SG_DateTime); %Calculate month
SG_T.Day=day(SG_T.SG_DateTime); %Calcualte day
SG_T.Hour=hour(SG_T.SG_DateTime); %Calculate hour
SG_T.Datenum=datenum(dateshift(SG_T.SG_DateTime,'start','day')); %Calculate day as number
SG_T.Season=SG_T.Month==11 | SG_T.Month==12 | SG_T.Month==1 | SG_T.Month==2 | SG_T.Month==3 | SG_T.Month==4; %So season =1 if between Nov and April (wet season), will equal 0 if the dry season (May to Oct)

%Now apply grpstats
SG_AH=grpstats(SG_T(:,2:end),'Hour','mean'); %Calculates hourly average for all (NOTE THIS MIGHT NOT ALWAYS BE SUITABLE) Also note 2:end needed as cannot average DateTime
SG_A_season=grpstats(SG_T(:,2:end),'Season','mean'); %Average for each season
SG_A_all=mean(SG_T{:,2:end}); %Mean of all columns
SG_S_all=sum(SG_T{:,2:end}); %Sum of all columns

%Matrix - differences are change - standard
SG_Smelt_A_all_grid=reshape(SG_A_all(1:35),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
SG_Smelt_A_all_grid_diff=SG_Smelt_A_all_grid-SG_Smelt_A_all_grid(1,3);
SG_Imelt_A_all_grid=reshape(SG_A_all(36:70),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
SG_Imelt_A_all_grid_diff=SG_Imelt_A_all_grid-SG_Imelt_A_all_grid(1,3);
SG_Tmelt_A_all_grid=reshape(SG_A_all(71:105),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
SG_Tmelt_A_all_grid_diff=SG_Tmelt_A_all_grid-SG_Tmelt_A_all_grid(1,3);
SG_Sub_A_all_grid=reshape(SG_A_all(176:210),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
SG_Sub_A_all_grid_diff=SG_Sub_A_all_grid-SG_Sub_A_all_grid(1,3);
SG_Abl_A_all_grid=reshape(SG_A_all(211:245),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
SG_Abl_A_all_grid_diff=SG_Abl_A_all_grid-SG_Abl_A_all_grid(1,3);

SG_Pr_sno_S_all_grid=reshape(SG_S_all(351:385),[7,5]); %This is Ta x Pr change At(1,3) is the standard run. This is the sum of Pr_sno
SG_Pr_liq_S_all_grid=reshape(SG_S_all(316:350),[7,5]); %This is Ta x Pr change At(1,3) is the standard run. This is the sum of Pr_liq
SG_percsno_grid=(SG_Pr_sno_S_all_grid./(SG_Pr_sno_S_all_grid+SG_Pr_liq_S_all_grid))*100; %Percentage snow of total precipitation

%Difference in mass balance
SG_D_all=SG_T{end,2:end}-SG_T{1,2:end}; %Difference between end and beginning for mass balances
SG_ICE_D_all_grid=reshape(SG_D_all(386:420),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
SG_ICE_D_all_grid_diff=SG_ICE_D_all_grid-SG_ICE_D_all_grid(1,3);
SG_ICE_D_all_grid_diff_m=SG_ICE_D_all_grid_diff./1000; %Now in m w.e.
SG_SWE_D_all_grid=reshape(SG_D_all(421:455),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
SG_SWE_D_all_grid_diff=SG_SWE_D_all_grid-SG_SWE_D_all_grid(1,3);
SG_SWE_D_all_grid_diff_m=SG_SWE_D_all_grid_diff./1000; %Now in m w.e.
SG_ISWE_D_all_grid=reshape(SG_D_all(456:490),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
SG_ISWE_D_all_grid_diff=SG_ISWE_D_all_grid-SG_ISWE_D_all_grid(1,3);
SG_ISWE_D_all_grid_diff_m=SG_ISWE_D_all_grid_diff./1000; %Now in m w.e.
SG_ISWE_D_all_grid_diff_myr=SG_ISWE_D_all_grid_diff_m./(SG_datasize/(365*24)); %So a rough m/yr based on the number of hours in the record

%=======================================================
%Artesonraju
%=======================================================
%Useful things
AG_datasize=size(AG(1).run.Date,1);

%Create matrix for each variable
AG_vars = [AG(1:35).run];
AG_Smelt=[AG_vars(1:35).Smelt]; %yes! this pulls out both Smelts into a table
AG_Imelt=[AG_vars(1:35).Imelt];
AG_Tmelt=AG_Smelt+AG_Imelt; %Good this works - total melt
AG_ESN=[AG_vars(1:35).ESN];
AG_EICE=[AG_vars(1:35).EICE];
AG_Sub=AG_ESN+AG_EICE; %Sublimation
AG_Abl=AG_Smelt+AG_Imelt+AG_ESN+AG_EICE; %Ablation
AG_Ta=[AG_vars(1:35).Ta];
AG_Pr=[AG_vars(1:35).Pr];
AG_Pr_liq=[AG_vars(1:35).Pr_liq];
AG_Pr_sno=[AG_vars(1:35).Pr_sno];
AG_ICE=[AG_vars(1:35).ICE];
AG_SWE=[AG_vars(1:35).SWE];
AG_ISWE=AG_ICE+AG_SWE;

AG_T=[AG_Smelt AG_Imelt AG_Tmelt AG_ESN AG_EICE AG_Sub AG_Abl AG_Ta AG_Pr AG_Pr_liq AG_Pr_sno AG_ICE AG_SWE AG_ISWE]; %Combine data into one matrix

%Create big table
AG_T=array2table(AG_T,'VariableNames',Table_names_long); %Note use of array2table when turning a matrix into a table
AG_DateTime=datetime(AG_vars(1).Date,'ConvertFrom','datenum');
AG_T = addvars(AG_T,AG_DateTime,'Before',1);
AG_T.Year=year(AG_T.AG_DateTime); %Calculate year
AG_T.Month=month(AG_T.AG_DateTime); %Calculate month
AG_T.Day=day(AG_T.AG_DateTime); %Calcualte day
AG_T.Hour=hour(AG_T.AG_DateTime); %Calculate hour
AG_T.Datenum=datenum(dateshift(AG_T.AG_DateTime,'start','day')); %Calculate day as number
AG_T.Season=AG_T.Month==11 | AG_T.Month==12 | AG_T.Month==1 | AG_T.Month==2 | AG_T.Month==3 | AG_T.Month==4; %So season =1 if between Nov and April (wet season), will equal 0 if the dry season (May to Oct)

%Now apply grpstats
AG_AH=grpstats(AG_T(:,2:end),'Hour','mean'); %Calculates hourly average for all (NOTE THIS MIGHT NOT ALWAYS BE SUITABLE) Also note 2:end needed as cannot average DateTime
AG_A_season=grpstats(AG_T(:,2:end),'Season','mean'); %Average for each season
AG_A_all=mean(AG_T{:,2:end}); %Mean of all columns
AG_S_all=sum(AG_T{:,2:end}); %Sum of all columns

%Matrix - differences are change - standard
AG_Smelt_A_all_grid=reshape(AG_A_all(1:35),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
AG_Smelt_A_all_grid_diff=AG_Smelt_A_all_grid-AG_Smelt_A_all_grid(1,3);
AG_Imelt_A_all_grid=reshape(AG_A_all(36:70),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
AG_Imelt_A_all_grid_diff=AG_Imelt_A_all_grid-AG_Imelt_A_all_grid(1,3);
AG_Tmelt_A_all_grid=reshape(AG_A_all(71:105),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
AG_Tmelt_A_all_grid_diff=AG_Tmelt_A_all_grid-AG_Tmelt_A_all_grid(1,3);
AG_Sub_A_all_grid=reshape(AG_A_all(176:210),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
AG_Sub_A_all_grid_diff=AG_Sub_A_all_grid-AG_Sub_A_all_grid(1,3);
AG_Abl_A_all_grid=reshape(AG_A_all(211:245),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
AG_Abl_A_all_grid_diff=AG_Abl_A_all_grid-AG_Abl_A_all_grid(1,3);

AG_Pr_sno_S_all_grid=reshape(AG_S_all(351:385),[7,5]); %This is Ta x Pr change At(1,3) is the standard run. This is the sum of Pr_sno
AG_Pr_liq_S_all_grid=reshape(AG_S_all(316:350),[7,5]); %This is Ta x Pr change At(1,3) is the standard run. This is the sum of Pr_liq
AG_percsno_grid=(AG_Pr_sno_S_all_grid./(AG_Pr_sno_S_all_grid+AG_Pr_liq_S_all_grid))*100; %Percentage snow of total precipitation

%Difference in mass balance
AG_D_all=AG_T{end,2:end}-AG_T{1,2:end}; %Difference between end and beginning for mass balances
AG_ICE_D_all_grid=reshape(AG_D_all(386:420),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
AG_ICE_D_all_grid_diff=AG_ICE_D_all_grid-AG_ICE_D_all_grid(1,3);
AG_ICE_D_all_grid_diff_m=AG_ICE_D_all_grid_diff./1000; %Now in m w.e.
AG_SWE_D_all_grid=reshape(AG_D_all(421:455),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
AG_SWE_D_all_grid_diff=AG_SWE_D_all_grid-AG_SWE_D_all_grid(1,3);
AG_SWE_D_all_grid_diff_m=AG_SWE_D_all_grid_diff./1000; %Now in m w.e.
AG_ISWE_D_all_grid=reshape(AG_D_all(456:490),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
AG_ISWE_D_all_grid_diff=AG_ISWE_D_all_grid-AG_ISWE_D_all_grid(1,3);
AG_ISWE_D_all_grid_diff_m=AG_ISWE_D_all_grid_diff./1000; %Now in m w.e.
AG_ISWE_D_all_grid_diff_myr=AG_ISWE_D_all_grid_diff_m./(AG_datasize/(365*24)); %So a rough m/yr based on the number of hours in the record

%======================================================
%Cuchillacocha
%======================================================
%Useful things
CG_datasize=size(CG(1).run.Date,1);

%Create matrix for each variable
CG_vars = [CG(1:35).run];
CG_Smelt=[CG_vars(1:35).Smelt]; %yes! this pulls out both Smelts into a table
CG_Imelt=[CG_vars(1:35).Imelt];
CG_Tmelt=CG_Smelt+CG_Imelt; %Good this works - total melt
CG_ESN=[CG_vars(1:35).ESN];
CG_EICE=[CG_vars(1:35).EICE];
CG_Sub=CG_ESN+CG_EICE; %Sublimation
CG_Abl=CG_Smelt+CG_Imelt+CG_ESN+CG_EICE; %Ablation
CG_Ta=[CG_vars(1:35).Ta];
CG_Pr=[CG_vars(1:35).Pr];
CG_Pr_liq=[CG_vars(1:35).Pr_liq];
CG_Pr_sno=[CG_vars(1:35).Pr_sno];
CG_ICE=[CG_vars(1:35).ICE];
CG_SWE=[CG_vars(1:35).SWE];
CG_ISWE=CG_ICE+CG_SWE;

CG_T=[CG_Smelt CG_Imelt CG_Tmelt CG_ESN CG_EICE CG_Sub CG_Abl CG_Ta CG_Pr CG_Pr_liq CG_Pr_sno CG_ICE CG_SWE CG_ISWE]; %Combine data into one matrix

%Create big table
CG_T=array2table(CG_T,'VariableNames',Table_names_long); %Note use of array2table when turning a matrix into a table
CG_DateTime=datetime(CG_vars(1).Date,'ConvertFrom','datenum');
CG_T = addvars(CG_T,CG_DateTime,'Before',1);
CG_T.Year=year(CG_T.CG_DateTime); %Calculate year
CG_T.Month=month(CG_T.CG_DateTime); %Calculate month
CG_T.Day=day(CG_T.CG_DateTime); %Calcualte day
CG_T.Hour=hour(CG_T.CG_DateTime); %Calculate hour
CG_T.Datenum=datenum(dateshift(CG_T.CG_DateTime,'start','day')); %Calculate day as number
CG_T.Season=CG_T.Month==11 | CG_T.Month==12 | CG_T.Month==1 | CG_T.Month==2 | CG_T.Month==3 | CG_T.Month==4; %So season =1 if between Nov and April (wet season), will equal 0 if the dry season (May to Oct)

%Now apply grpstats
CG_AH=grpstats(CG_T(:,2:end),'Hour','mean'); %Calculates hourly average for all (NOTE THIS MIGHT NOT ALWAYS BE SUITABLE) Also note 2:end needed as cannot average DateTime
CG_A_season=grpstats(CG_T(:,2:end),'Season','mean'); %Average for each season
CG_A_all=mean(CG_T{:,2:end}); %Mean of all columns
CG_S_all=sum(CG_T{:,2:end}); %Sum of all columns

%Matrix - differences are change - standard
CG_Smelt_A_all_grid=reshape(CG_A_all(1:35),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
CG_Smelt_A_all_grid_diff=CG_Smelt_A_all_grid-CG_Smelt_A_all_grid(1,3);
CG_Imelt_A_all_grid=reshape(CG_A_all(36:70),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
CG_Imelt_A_all_grid_diff=CG_Imelt_A_all_grid-CG_Imelt_A_all_grid(1,3);
CG_Tmelt_A_all_grid=reshape(CG_A_all(71:105),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
CG_Tmelt_A_all_grid_diff=CG_Tmelt_A_all_grid-CG_Tmelt_A_all_grid(1,3);
CG_Sub_A_all_grid=reshape(CG_A_all(176:210),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
CG_Sub_A_all_grid_diff=CG_Sub_A_all_grid-CG_Sub_A_all_grid(1,3);
CG_Abl_A_all_grid=reshape(CG_A_all(211:245),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
CG_Abl_A_all_grid_diff=CG_Abl_A_all_grid-CG_Abl_A_all_grid(1,3);

CG_Pr_sno_S_all_grid=reshape(CG_S_all(351:385),[7,5]); %This is Ta x Pr change At(1,3) is the standard run. This is the sum of Pr_sno
CG_Pr_liq_S_all_grid=reshape(CG_S_all(316:350),[7,5]); %This is Ta x Pr change At(1,3) is the standard run. This is the sum of Pr_liq
CG_percsno_grid=(CG_Pr_sno_S_all_grid./(CG_Pr_sno_S_all_grid+CG_Pr_liq_S_all_grid))*100; %Percentage snow of total precipitation

%Difference in mass balance
CG_D_all=CG_T{end,2:end}-CG_T{1,2:end}; %Difference between end and beginning for mass balances
CG_ICE_D_all_grid=reshape(CG_D_all(386:420),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
CG_ICE_D_all_grid_diff=CG_ICE_D_all_grid-CG_ICE_D_all_grid(1,3);
CG_ICE_D_all_grid_diff_m=CG_ICE_D_all_grid_diff./1000; %Now in m w.e.
CG_SWE_D_all_grid=reshape(CG_D_all(421:455),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
CG_SWE_D_all_grid_diff=CG_SWE_D_all_grid-CG_SWE_D_all_grid(1,3);
CG_SWE_D_all_grid_diff_m=CG_SWE_D_all_grid_diff./1000; %Now in m w.e.
CG_ISWE_D_all_grid=reshape(CG_D_all(456:490),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
CG_ISWE_D_all_grid_diff=CG_ISWE_D_all_grid-CG_ISWE_D_all_grid(1,3);
CG_ISWE_D_all_grid_diff_m=CG_ISWE_D_all_grid_diff./1000; %Now in m w.e.
CG_ISWE_D_all_grid_diff_myr=CG_ISWE_D_all_grid_diff_m./(CG_datasize/(365*24)); %So a rough m/yr based on the number of hours in the record

%=======================================================
%Quisoquipina
%=======================================================

%Useful things
QQG_datasize=size(QQG(1).run.Date,1);

%Create matrix for each variable
QQG_vars = [QQG(1:35).run];
QQG_Smelt=[QQG_vars(1:35).Smelt]; %yes! this pulls out both Smelts into a table
QQG_Imelt=[QQG_vars(1:35).Imelt];
QQG_Tmelt=QQG_Smelt+QQG_Imelt; %Good this works - total melt
QQG_ESN=[QQG_vars(1:35).ESN];
QQG_EICE=[QQG_vars(1:35).EICE];
QQG_Sub=QQG_ESN+QQG_EICE; %Sublimation
QQG_Abl=QQG_Smelt+QQG_Imelt+QQG_ESN+QQG_EICE; %Ablation
QQG_Ta=[QQG_vars(1:35).Ta];
QQG_Pr=[QQG_vars(1:35).Pr];
QQG_Pr_liq=[QQG_vars(1:35).Pr_liq];
QQG_Pr_sno=[QQG_vars(1:35).Pr_sno];
QQG_ICE=[QQG_vars(1:35).ICE];
QQG_SWE=[QQG_vars(1:35).SWE];
QQG_ISWE=QQG_ICE+QQG_SWE;

QQG_T=[QQG_Smelt QQG_Imelt QQG_Tmelt QQG_ESN QQG_EICE QQG_Sub QQG_Abl QQG_Ta QQG_Pr QQG_Pr_liq QQG_Pr_sno QQG_ICE QQG_SWE QQG_ISWE]; %Combine data into one matrix

%Create big table
QQG_T=array2table(QQG_T,'VariableNames',Table_names_long); %Note use of array2table when turning a matrix into a table
QQG_DateTime=datetime(QQG_vars(1).Date,'ConvertFrom','datenum');
QQG_T = addvars(QQG_T,QQG_DateTime,'Before',1);
QQG_T.Year=year(QQG_T.QQG_DateTime); %Calculate year
QQG_T.Month=month(QQG_T.QQG_DateTime); %Calculate month
QQG_T.Day=day(QQG_T.QQG_DateTime); %Calcualte day
QQG_T.Hour=hour(QQG_T.QQG_DateTime); %Calculate hour
QQG_T.Datenum=datenum(dateshift(QQG_T.QQG_DateTime,'start','day')); %Calculate day as number
QQG_T.Season=QQG_T.Month==11 | QQG_T.Month==12 | QQG_T.Month==1 | QQG_T.Month==2 | QQG_T.Month==3 | QQG_T.Month==4; %So season =1 if between Nov and April (wet season), will equal 0 if the dry season (May to Oct)

%Now apply grpstats
QQG_AH=grpstats(QQG_T(:,2:end),'Hour','mean'); %Calculates hourly average for all (NOTE THIS MIGHT NOT ALWAYS BE SUITABLE) Also note 2:end needed as cannot average DateTime
QQG_A_season=grpstats(QQG_T(:,2:end),'Season','mean'); %Average for each season
QQG_A_all=mean(QQG_T{:,2:end}); %Mean of all columns
QQG_S_all=sum(QQG_T{:,2:end}); %Sum of all columns

%Matrix - differences are change - standard
QQG_Smelt_A_all_grid=reshape(QQG_A_all(1:35),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
QQG_Smelt_A_all_grid_diff=QQG_Smelt_A_all_grid-QQG_Smelt_A_all_grid(1,3);
QQG_Imelt_A_all_grid=reshape(QQG_A_all(36:70),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
QQG_Imelt_A_all_grid_diff=QQG_Imelt_A_all_grid-QQG_Imelt_A_all_grid(1,3);
QQG_Tmelt_A_all_grid=reshape(QQG_A_all(71:105),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
QQG_Tmelt_A_all_grid_diff=QQG_Tmelt_A_all_grid-QQG_Tmelt_A_all_grid(1,3);
QQG_Sub_A_all_grid=reshape(QQG_A_all(176:210),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
QQG_Sub_A_all_grid_diff=QQG_Sub_A_all_grid-QQG_Sub_A_all_grid(1,3);
QQG_Abl_A_all_grid=reshape(QQG_A_all(211:245),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
QQG_Abl_A_all_grid_diff=QQG_Abl_A_all_grid-QQG_Abl_A_all_grid(1,3);

QQG_Pr_sno_S_all_grid=reshape(QQG_S_all(351:385),[7,5]); %This is Ta x Pr change At(1,3) is the standard run. This is the sum of Pr_sno
QQG_Pr_liq_S_all_grid=reshape(QQG_S_all(316:350),[7,5]); %This is Ta x Pr change At(1,3) is the standard run. This is the sum of Pr_liq
QQG_percsno_grid=(QQG_Pr_sno_S_all_grid./(QQG_Pr_sno_S_all_grid+QQG_Pr_liq_S_all_grid))*100; %Percentage snow of total precipitation

%Difference in mass balance
QQG_D_all=QQG_T{end,2:end}-QQG_T{1,2:end}; %Difference between end and beginning for mass balances
QQG_ICE_D_all_grid=reshape(QQG_D_all(386:420),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
QQG_ICE_D_all_grid_diff=QQG_ICE_D_all_grid-QQG_ICE_D_all_grid(1,3);
QQG_ICE_D_all_grid_diff_m=QQG_ICE_D_all_grid_diff./1000; %Now in m w.e.
QQG_SWE_D_all_grid=reshape(QQG_D_all(421:455),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
QQG_SWE_D_all_grid_diff=QQG_SWE_D_all_grid-QQG_SWE_D_all_grid(1,3);
QQG_SWE_D_all_grid_diff_m=QQG_SWE_D_all_grid_diff./1000; %Now in m w.e.
QQG_ISWE_D_all_grid=reshape(QQG_D_all(456:490),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
QQG_ISWE_D_all_grid_diff=QQG_ISWE_D_all_grid-QQG_ISWE_D_all_grid(1,3);
QQG_ISWE_D_all_grid_diff_m=QQG_ISWE_D_all_grid_diff./1000; %Now in m w.e.
QQG_ISWE_D_all_grid_diff_myr=QQG_ISWE_D_all_grid_diff_m./(QQG_datasize/(365*24)); %So a rough m/yr based on the number of hours in the record

%======================================================
%Quelccaya
%======================================================
%Useful things
QIC_datasize=size(QIC(1).run.Date,1);

%Create matrix for each variable
QIC_vars = [QIC(1:35).run];
QIC_Smelt=[QIC_vars(1:35).Smelt]; %yes! this pulls out both Smelts into a table
QIC_Imelt=[QIC_vars(1:35).Imelt];
QIC_Tmelt=QIC_Smelt+QIC_Imelt; %Good this works - total melt
QIC_ESN=[QIC_vars(1:35).ESN];
QIC_EICE=[QIC_vars(1:35).EICE];
QIC_Sub=QIC_ESN+QIC_EICE; %Sublimation
QIC_Abl=QIC_Smelt+QIC_Imelt+QIC_ESN+QIC_EICE; %Ablation
QIC_Ta=[QIC_vars(1:35).Ta];
QIC_Pr=[QIC_vars(1:35).Pr];
QIC_Pr_liq=[QIC_vars(1:35).Pr_liq];
QIC_Pr_sno=[QIC_vars(1:35).Pr_sno];
QIC_ICE=[QIC_vars(1:35).ICE];
QIC_SWE=[QIC_vars(1:35).SWE];
QIC_ISWE=QIC_ICE+QIC_SWE;

QIC_T=[QIC_Smelt QIC_Imelt QIC_Tmelt QIC_ESN QIC_EICE QIC_Sub QIC_Abl QIC_Ta QIC_Pr QIC_Pr_liq QIC_Pr_sno QIC_ICE QIC_SWE QIC_ISWE]; %Combine data into one matrix

%Create big table
QIC_T=array2table(QIC_T,'VariableNames',Table_names_long); %Note use of array2table when turning a matrix into a table
QIC_DateTime=datetime(QIC_vars(1).Date,'ConvertFrom','datenum');
QIC_T = addvars(QIC_T,QIC_DateTime,'Before',1);
QIC_T.Year=year(QIC_T.QIC_DateTime); %Calculate year
QIC_T.Month=month(QIC_T.QIC_DateTime); %Calculate month
QIC_T.Day=day(QIC_T.QIC_DateTime); %Calcualte day
QIC_T.Hour=hour(QIC_T.QIC_DateTime); %Calculate hour
QIC_T.Datenum=datenum(dateshift(QIC_T.QIC_DateTime,'start','day')); %Calculate day as number
QIC_T.Season=QIC_T.Month==11 | QIC_T.Month==12 | QIC_T.Month==1 | QIC_T.Month==2 | QIC_T.Month==3 | QIC_T.Month==4; %So season =1 if between Nov and April (wet season), will equal 0 if the dry season (May to Oct)

%Now apply grpstats
QIC_AH=grpstats(QIC_T(:,2:end),'Hour','mean'); %Calculates hourly average for all (NOTE THIS MIGHT NOT ALWAYS BE SUITABLE) Also note 2:end needed as cannot average DateTime
QIC_A_season=grpstats(QIC_T(:,2:end),'Season','mean'); %Average for each season
QIC_A_all=mean(QIC_T{:,2:end}); %Mean of all columns
QIC_S_all=sum(QIC_T{:,2:end}); %Sum of all columns

%Matrix - differences are change - standard
QIC_Smelt_A_all_grid=reshape(QIC_A_all(1:35),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
QIC_Smelt_A_all_grid_diff=QIC_Smelt_A_all_grid-QIC_Smelt_A_all_grid(1,3);
QIC_Imelt_A_all_grid=reshape(QIC_A_all(36:70),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
QIC_Imelt_A_all_grid_diff=QIC_Imelt_A_all_grid-QIC_Imelt_A_all_grid(1,3);
QIC_Tmelt_A_all_grid=reshape(QIC_A_all(71:105),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
QIC_Tmelt_A_all_grid_diff=QIC_Tmelt_A_all_grid-QIC_Tmelt_A_all_grid(1,3);
QIC_Sub_A_all_grid=reshape(QIC_A_all(176:210),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
QIC_Sub_A_all_grid_diff=QIC_Sub_A_all_grid-QIC_Sub_A_all_grid(1,3);
QIC_Abl_A_all_grid=reshape(QIC_A_all(211:245),[7,5]); %This is Ta x Pr change At(1,3) is the standard run
QIC_Abl_A_all_grid_diff=QIC_Abl_A_all_grid-QIC_Abl_A_all_grid(1,3);

QIC_Pr_sno_S_all_grid=reshape(QIC_S_all(351:385),[7,5]); %This is Ta x Pr change At(1,3) is the standard run. This is the sum of Pr_sno
QIC_Pr_liq_S_all_grid=reshape(QIC_S_all(316:350),[7,5]); %This is Ta x Pr change At(1,3) is the standard run. This is the sum of Pr_liq
QIC_percsno_grid=(QIC_Pr_sno_S_all_grid./(QIC_Pr_sno_S_all_grid+QIC_Pr_liq_S_all_grid))*100; %Percentage snow of total precipitation

%Difference in mass balance
QIC_D_all=QIC_T{end,2:end}-QIC_T{1,2:end}; %Difference between end and beginning for mass balances
QIC_ICE_D_all_grid=reshape(QIC_D_all(386:420),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
QIC_ICE_D_all_grid_diff=QIC_ICE_D_all_grid-QIC_ICE_D_all_grid(1,3);
QIC_ICE_D_all_grid_diff_m=QIC_ICE_D_all_grid_diff./1000; %Now in m w.e.
QIC_SWE_D_all_grid=reshape(QIC_D_all(421:455),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
QIC_SWE_D_all_grid_diff=QIC_SWE_D_all_grid-QIC_SWE_D_all_grid(1,3);
QIC_SWE_D_all_grid_diff_m=QIC_SWE_D_all_grid_diff./1000; %Now in m w.e.
QIC_ISWE_D_all_grid=reshape(QIC_D_all(456:490),[7,5]);%This is Ta x Pr change At(1,3) is the standard run
QIC_ISWE_D_all_grid_diff=QIC_ISWE_D_all_grid-QIC_ISWE_D_all_grid(1,3);
QIC_ISWE_D_all_grid_diff_m=QIC_ISWE_D_all_grid_diff./1000; %Now in m w.e.
QIC_ISWE_D_all_grid_diff_myr=QIC_ISWE_D_all_grid_diff_m./(QIC_datasize/(365*24)); %So a rough m/yr based on the number of hours in the record

%=======================================================
%Figures
%========================================================

%========== Plot with just varying Ta and Pr ===========================

Nice_colours=[0.502,0,0;0,0,0.502;0.961,0.510,0.188;0.902,0.098,0.294;0,0.510,0.784;0.863,0.745,1;0.980,0.745,0.831;0.275,0.941,0.941];
%For info ax=subplot(8,4,1);  ax.ColorOrder = Nice_colours; hold
Fluxes_colours=[Nice_colours(2,:); Nice_colours(4,:); Nice_colours(3,:); Nice_colours(5,:); Nice_colours(6,:); Nice_colours(7,:); Nice_colours(8,:)];
All_names=["SG","AG","CG","QQG","QIC"];

figure2=figure();
t = tiledlayout(2,5,'TileSpacing','Compact');
ax=nexttile(1); ax.ColorOrder = Nice_colours; hold on;
plot(Ta_change_short,SG_Smelt_A_all_grid_diff(:,3),'-d',Ta_change_short,AG_Smelt_A_all_grid_diff(:,3),'-d',Ta_change_short,CG_Smelt_A_all_grid_diff(:,3),'-d',Ta_change_short,QQG_Smelt_A_all_grid_diff(:,3),'-d',Ta_change_short,QIC_Smelt_A_all_grid_diff(:,3),'-d');
title('Snow melt'); ylabel('Snow melt diff. (mm w.e. h^-^1)'); xlabel('Air temperature diff. (\circC)'); 

ax=nexttile(2); ax.ColorOrder = Nice_colours; hold on;
plot(Ta_change_short,SG_Imelt_A_all_grid_diff(:,3),'-d',Ta_change_short,AG_Imelt_A_all_grid_diff(:,3),'-d',Ta_change_short,CG_Imelt_A_all_grid_diff(:,3),'-d',Ta_change_short,QQG_Imelt_A_all_grid_diff(:,3),'-d',Ta_change_short,QIC_Imelt_A_all_grid_diff(:,3),'-d');
title('Ice melt'); ylabel('Ice melt diff. (mm w.e. h^-^1)'); xlabel('Air temperature diff. (\circC)');

ax=nexttile(3); ax.ColorOrder = Nice_colours; hold on;
plot(Ta_change_short,SG_Tmelt_A_all_grid_diff(:,3),'-d',Ta_change_short,AG_Tmelt_A_all_grid_diff(:,3),'-d',Ta_change_short,CG_Tmelt_A_all_grid_diff(:,3),'-d',Ta_change_short,QQG_Tmelt_A_all_grid_diff(:,3),'-d',Ta_change_short,QIC_Tmelt_A_all_grid_diff(:,3),'-d');
title('Total melt'); ylabel('Total melt diff. (mm w.e. h^-^1)'); xlabel('Air temperature diff. (\circC)');

ax=nexttile(4); ax.ColorOrder = Nice_colours; hold on;
plot(Ta_change_short,SG_Sub_A_all_grid_diff(:,3),'-d',Ta_change_short,AG_Sub_A_all_grid_diff(:,3),'-d',Ta_change_short,CG_Sub_A_all_grid_diff(:,3),'-d',Ta_change_short,QQG_Sub_A_all_grid_diff(:,3),'-d',Ta_change_short,QIC_Sub_A_all_grid_diff(:,3),'-d');
title('Sublimation'); ylabel('Sublimation diff. (mm w.e. h^-^1)'); xlabel('Air temperature diff. (\circC)');

ax=nexttile(5); ax.ColorOrder = Nice_colours; hold on;
plot(Ta_change_short,SG_Abl_A_all_grid_diff(:,3),'-d',Ta_change_short,AG_Abl_A_all_grid_diff(:,3),'-d',Ta_change_short,CG_Abl_A_all_grid_diff(:,3),'-d',Ta_change_short,QQG_Abl_A_all_grid_diff(:,3),'-d',Ta_change_short,QIC_Abl_A_all_grid_diff(:,3),'-d');
title('Ablation'); ylabel('Ablation diff. (mm w.e. h^-^1)'); xlabel('Air temperature diff. (\circC)');

ax=nexttile(6); ax.ColorOrder = Nice_colours; hold on;
plot(Pr_change_name_short,SG_Smelt_A_all_grid_diff(1,:),'-d',Pr_change_name_short,AG_Smelt_A_all_grid_diff(1,:),'-d',Pr_change_name_short,CG_Smelt_A_all_grid_diff(1,:),'-d',Pr_change_name_short,QQG_Smelt_A_all_grid_diff(1,:),'-d',Pr_change_name_short,QIC_Smelt_A_all_grid_diff(1,:),'-d');
title('Snow melt'); ylabel('Snow melt diff. (mm w.e. h^-^1)'); xlabel('Precipitation diff. (%)'); xticks(Pr_change_name_short);

ax=nexttile(7); ax.ColorOrder = Nice_colours; hold on;
plot(Pr_change_name_short,SG_Imelt_A_all_grid_diff(1,:),'-d',Pr_change_name_short,AG_Imelt_A_all_grid_diff(1,:),'-d',Pr_change_name_short,CG_Imelt_A_all_grid_diff(1,:),'-d',Pr_change_name_short,QQG_Imelt_A_all_grid_diff(1,:),'-d',Pr_change_name_short,QIC_Imelt_A_all_grid_diff(1,:),'-d');
title('Ice melt'); ylabel('Ice melt diff. (mm w.e. h^-^1)'); xlabel('Precipitation diff. (%)');xticks(Pr_change_name_short);

ax=nexttile(8); ax.ColorOrder = Nice_colours; hold on;
plot(Pr_change_name_short,SG_Tmelt_A_all_grid_diff(1,:),'-d',Pr_change_name_short,AG_Tmelt_A_all_grid_diff(1,:),'-d',Pr_change_name_short,CG_Tmelt_A_all_grid_diff(1,:),'-d',Pr_change_name_short,QQG_Tmelt_A_all_grid_diff(1,:),'-d',Pr_change_name_short,QIC_Tmelt_A_all_grid_diff(1,:),'-d');
title('Total melt'); ylabel('Total melt diff. (mm w.e. h^-^1)'); xlabel('Precipitation diff. (%)');xticks(Pr_change_name_short);
legend(All_names,'Location','SouthOutside','Orientation','Horizontal');

ax=nexttile(9); ax.ColorOrder = Nice_colours; hold on;
plot(Pr_change_name_short,SG_Sub_A_all_grid_diff(1,:),'-d',Pr_change_name_short,AG_Sub_A_all_grid_diff(1,:),'-d',Pr_change_name_short,CG_Sub_A_all_grid_diff(1,:),'-d',Pr_change_name_short,QQG_Sub_A_all_grid_diff(1,:),'-d',Pr_change_name_short,QIC_Sub_A_all_grid_diff(1,:),'-d');
title('Sublimation'); ylabel('Sublimation diff. (mm w.e. h^-^1)'); xlabel('Precipitation diff. (%)');xticks(Pr_change_name_short);

ax=nexttile(10); ax.ColorOrder = Nice_colours; hold on;
plot(Pr_change_name_short,SG_Abl_A_all_grid_diff(1,:),'-d',Pr_change_name_short,AG_Abl_A_all_grid_diff(1,:),'-d',Pr_change_name_short,CG_Abl_A_all_grid_diff(1,:),'-d',Pr_change_name_short,QQG_Abl_A_all_grid_diff(1,:),'-d',Pr_change_name_short,QIC_Abl_A_all_grid_diff(1,:),'-d');
title('Ablation'); ylabel('Ablation diff. (mm w.e. h^-^1)'); xlabel('Precipitation diff. (%)');xticks(Pr_change_name_short);

% Create textbox
annotation(figure2,'textbox',[0.0590399619410085 0.934905660377357 0.0185052331113225 0.0499999999999987],'String',{'a'},'LineStyle','none','FontSize',14,'FitBoxToText','off');
annotation(figure2,'textbox',[0.056185537583254 0.481132075471696 0.0185052331113224 0.0499999999999987],'String','b','LineStyle','none','FontSize',14,'FitBoxToText','off');


%============= Mass change heat map ===============================
load('ISWE_color2.mat','ISWE_map'); %use bluewhitered_heatmap

figure3=figure(); 
t = tiledlayout(2,6,'TileSpacing','Compact');
nexttile(1,[1 2]);
h1=heatmap(Pr_change_name_short,Ta_change_short,SG_ISWE_D_all_grid_diff_myr);
h1.CellLabelFormat = '%.1f'; %Shows numbers to 1 decimal places
title('SG');
h1.ColorLimits = [-14 2];
h1.Colormap = ISWE_map; colorbar;
nexttile(3,[1 2]);
h1=heatmap(Pr_change_name_short,Ta_change_short,AG_ISWE_D_all_grid_diff_myr);
h1.CellLabelFormat = '%.1f'; %Shows numbers to 1 decimal places
title('AG');
h1.ColorLimits = [-14 2];
h1.Colormap = ISWE_map; colorbar; 
nexttile(5, [1 2]);
h1=heatmap(Pr_change_name_short,Ta_change_short,CG_ISWE_D_all_grid_diff_myr);
h1.CellLabelFormat = '%.1f'; %Shows numbers to 1 decimal places
title('CG');
h1.ColorLimits = [-14 2];
h1.Colormap = ISWE_map; colorbar; 
nexttile(8, [1 2]);
h1=heatmap(Pr_change_name_short,Ta_change_short,QQG_ISWE_D_all_grid_diff_myr);
h1.CellLabelFormat = '%.1f'; %Shows numbers to 1 decimal places
title('QQG');
h1.ColorLimits = [-14 2];
h1.Colormap = ISWE_map; colorbar; 
nexttile(10,[1 2]);
h1=heatmap(Pr_change_name_short,Ta_change_short,QIC_ISWE_D_all_grid_diff_myr);
h1.CellLabelFormat = '%.1f'; %Shows numbers to 1 decimal places
title('QIC');
h1.ColorLimits = [-14 2];
h1.Colormap = ISWE_map; colorbar; 
xlabel(t,'Precipitation difference (%)','FontSize',14); %Global titles
ylabel(t,'Air temperature difference (\circC)','FontSize',14);






