%Script to analyse all the T and C outputs for the Peruvian weather stations. 
%Should have run the validation code already. 

%By Cat Fyffe 22/09/2020 %Edited 22/08/2021 for input into EIDC

%Add path
addpath('C:\Users\catri\Dropbox\Northumbria_University\PeruGROWS\Energy_balance_paper\Data_for_paper\T_and_C_analysis\Main_runs_for_analysis');
%Paths for climate change data
addpath('C:\Users\catri\Dropbox\Northumbria_University\PeruGROWS\Modelling\T_and_C_HIMAL\Outputs_Peru\SG_outputs');
addpath('C:\Users\catri\Dropbox\Northumbria_University\PeruGROWS\Modelling\T_and_C_HIMAL\Outputs_Peru\AG_outputs');
addpath('C:\Users\catri\Dropbox\Northumbria_University\PeruGROWS\Modelling\T_and_C_HIMAL\Outputs_Peru\CG_outputs');
addpath('C:\Users\catri\Dropbox\Northumbria_University\PeruGROWS\Modelling\T_and_C_HIMAL\Outputs_Peru\QQG_outputs');
addpath('C:\Users\catri\Dropbox\Northumbria_University\PeruGROWS\Modelling\T_and_C_HIMAL\Outputs_Peru\QIC_outputs');
addpath('C:\Users\catri\Dropbox\Northumbria_University\PeruGROWS\Energy_balance_paper\Data_for_paper\T_and_C_analysis\Main_analysis');

%Bring in the data
SG=load('SG_Short_k.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Rsw_out','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ');
AG=load('AG_Run_ax.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Rsw_out','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ');
CG=load('CG_Run_z.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ');
QIC=load('QIC_Run_rosTa_510Pr_510.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ'); %This has the updated change to Snowpacks
QQG=load('QQG_Run_g_newval.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Rsw_out','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ');

% %Going to run under a Ta increase to compare
% SG=load('SG_Run_Ta_2Pr_0.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Rsw_out','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ');
% AG=load('AG_Run_Ta_2Pr_0.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Rsw_out','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ');
% CG=load('CG_Run_Ta_2Pr_0.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ');
% QIC=load('QIC_Run_Ta_2Pr_0.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ'); %This has the updated change to Snowpacks
% QQG=load('QQG_Run_Ta_2Pr_0.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Rsw_out','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ');
% %Note run validation code before this

% %Going to run with climate change standard run
% SG=load('SG_Run_Ta_0Pr_0.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Rsw_out','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ');
% AG=load('AG_Run_Ta_0Pr_0.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Rsw_out','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ');
% CG=load('CG_Run_Ta_0Pr_0.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ');
% QIC=load('QIC_Run_Ta_0Pr_0.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ'); %This has the updated change to Snowpacks
% QQG=load('QQG_Run_Ta_0Pr_0.mat','DateTime','ALB','Csno','ea','EICE','ESN','G','H','ICE','ICE_D','Imelt','IP_wc','L_out','Latm','N','Pr','Pr_liq','Pr_sno','Pre','QE','Qfm','Qv','Rn','ros','Rsw','Rsw_out','Smelt','SND','SR_50','SWE','Ta','Tice','Ts','U','Ws','dQ');
% %Note run validation code before this

%Helpful variables
SG_datasize=size(SG.DateTime,1);
AG_datasize=size(AG.DateTime,1);
CG_datasize=size(CG.DateTime,1);
QIC_datasize=size(QIC.DateTime,1);
QQG_datasize=size(QQG.DateTime,1);
SG_el=load('SG_Short_k.mat','Zbas');SG_el=SG_el.Zbas;
AG_el=load('AG_Run_ax.mat','Zbas');AG_el=AG_el.Zbas;
CG_el=load('CG_Run_z.mat','Zbas'); CG_el=CG_el.Zbas;
QIC_el=load('QIC_Run_rosTa_510Pr_510.mat','Zbas');QIC_el=QIC_el.Zbas;
QQG_el=load('QQG_Run_g_newval.mat','Zbas');QQG_el=QQG_el.Zbas;
All_el=[SG_el AG_el CG_el QQG_el QIC_el];

%Sort out radiation
%Rsw_out_mod - Outgoing Shortwave
SG.Rsw_out_mod=SG.Rsw.*SG.ALB;
SG.Rsw_out_mod(isnan(SG.Rsw_out_mod)==1)=0; %Replace NaN with 0
SG.Rsw_out_mod=-SG.Rsw_out_mod; %Make -ve as outgoing
AG.Rsw_out_mod=AG.Rsw.*AG.ALB;
AG.Rsw_out_mod(isnan(AG.Rsw_out_mod)==1)=0; %Replace NaN with 0
AG.Rsw_out_mod=-AG.Rsw_out_mod; %Make -ve as outgoing
CG.Rsw_out_mod=CG.Rsw.*CG.ALB;
CG.Rsw_out_mod(isnan(CG.Rsw_out_mod)==1)=0; %Replace NaN with 0
CG.Rsw_out_mod=-CG.Rsw_out_mod; %Make -ve as outgoing
QIC.Rsw_out_mod=QIC.Rsw.*QIC.ALB;
QIC.Rsw_out_mod(isnan(QIC.Rsw_out_mod)==1)=0; %Replace NaN with 0
QIC.Rsw_out_mod=-QIC.Rsw_out_mod; %Make -ve as outgoing
QQG.Rsw_out_mod=QQG.Rsw.*QQG.ALB;
QQG.Rsw_out_mod(isnan(QQG.Rsw_out_mod)==1)=0; %Replace NaN with 0
QQG.Rsw_out_mod=-QQG.Rsw_out_mod; %Make -ve as outgoing

%Calculate Incoming Longwave
SG.Latm_mod=zeros(SG_datasize,1);
AG.Latm_mod=zeros(AG_datasize,1);
CG.Latm_mod=zeros(CG_datasize,1);
QIC.Latm_mod=zeros(QIC_datasize,1);
QQG.Latm_mod=zeros(QQG_datasize,1);

for t=1:SG_datasize
    SG.Latm_mod(t)=Incoming_Longwave(SG.Ta(t),SG.ea(t),SG.N(t));
end
for t=1:AG_datasize
    AG.Latm_mod(t)=Incoming_Longwave(AG.Ta(t),AG.ea(t),AG.N(t));
end
for t=1:CG_datasize
    CG.Latm_mod(t)=Incoming_Longwave(CG.Ta(t),CG.ea(t),CG.N(t));
end
for t=1:QIC_datasize
    QIC.Latm_mod(t)=Incoming_Longwave(QIC.Ta(t),QIC.ea(t),QIC.N(t));
end
for t=1:QQG_datasize
    QQG.Latm_mod(t)=Incoming_Longwave(QQG.Ta(t),QQG.ea(t),QQG.N(t));
end

%Calculate Outgoing Longwave
a_sur=0.97;
sigmaSB = 5.6704e-8;
SG.Ts_k=SG.Ts+273.15;
SG.L_out_mod = a_sur.*sigmaSB.*(SG.Ts_k).^4; %Modelled Lout from modelled Ts
SG.L_out_mod=-SG.L_out_mod; %Make -ve as outgoing
AG.Ts_k=AG.Ts+273.15;
AG.L_out_mod = a_sur.*sigmaSB.*(AG.Ts_k).^4; %Modelled Lout from modelled Ts
AG.L_out_mod=-AG.L_out_mod; %Make -ve as outgoing
CG.Ts_k=CG.Ts+273.15;
CG.L_out_mod = a_sur.*sigmaSB.*(CG.Ts_k).^4; %Modelled Lout from modelled Ts
CG.L_out_mod=-CG.L_out_mod; %Make -ve as outgoing
QIC.Ts_k=QIC.Ts+273.15;
QIC.L_out_mod = a_sur.*sigmaSB.*(QIC.Ts_k).^4; %Modelled Lout from modelled Ts
QIC.L_out_mod=-QIC.L_out_mod; %Make -ve as outgoing
QQG.Ts_k=QQG.Ts+273.15;
QQG.L_out_mod = a_sur.*sigmaSB.*(QQG.Ts_k).^4; %Modelled Lout from modelled Ts
QQG.L_out_mod=-QQG.L_out_mod; %Make -ve as outgoing

%Net radiation
SG.SWnet=SG.Rsw+SG.Rsw_out_mod;
AG.SWnet=AG.Rsw+AG.Rsw_out_mod;
CG.SWnet=CG.Rsw+CG.Rsw_out_mod;
QIC.SWnet=QIC.Rsw+QIC.Rsw_out_mod;
QQG.SWnet=QQG.Rsw+QQG.Rsw_out_mod;
SG.LWnet=SG.Latm_mod+SG.L_out_mod;
AG.LWnet=AG.Latm_mod+AG.L_out_mod;
CG.LWnet=CG.Latm_mod+CG.L_out_mod;
QIC.LWnet=QIC.Latm_mod+QIC.L_out_mod;
QQG.LWnet=QQG.Latm_mod+QQG.L_out_mod;

%Change sign to convention of + towards surface and - away from surface
SG.H=-SG.H; AG.H=-AG.H; CG.H=-CG.H; QIC.H=-QIC.H; QQG.H=-QQG.H;%Sensible heat
SG.QE=-SG.QE; AG.QE=-AG.QE; CG.QE=-CG.QE; QIC.QE=-QIC.QE; QQG.QE=-QQG.QE; %Latent heat
SG.G=-SG.G; AG.G=-AG.G; CG.G=-CG.G; QIC.G=-QIC.G; QQG.G=-QQG.G; %ground heat flux
%Qfm, Qv and Rn are already positive


%Going to do the rest one at a time so easier to repeat

%For everyone
Lf= 1000*333.7; %% [J/Kg] %% Latent heat melt/freezing 
HyYear_threshold = datetime('01-11-0001','InputFormat','dd-MM-yyyy'); %Day to determine start of hydrological year

%====== Shallap Glacier ================

SG_T=struct2table(SG); %Create table
%Do some extra calculations
SG_T.Tmelt=SG_T.Smelt+SG_T.Imelt; %Total melt
SG_L= 1000.*(2501.3 - 2.361.*(SG_T.Ta)); %%% Latent heat vaporization/condensaition [J/kg]
SG_Ls= Lf+SG_L; %%% Latent Heat of sublimation  [J/kg]
SG_T.EICE_Wm2=(SG_T.EICE./3600).*SG_Ls;
SG_T.ESN_Wm2=(SG_T.ESN./3600).*SG_Ls;

%Sort time
SG_T.Year=year(SG_T.DateTime); %Calculate year
SG_T.Month=month(SG_T.DateTime); %Calculate month
SG_T.Day=day(SG_T.DateTime); %Calcualte day
SG_T.Hour=hour(SG_T.DateTime); %Calculate hour
SG_T.Datenum=datenum(dateshift(SG_T.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)
SG_Date=datetime(0001,SG_T.Month,SG_T.Day,SG_T.Hour,0,0);%This is the day and time but no year
SG_idhy=SG_Date>=HyYear_threshold; %So if after 1st Nov then = 1
SG_T.HyYear(SG_idhy)=SG_T.Year(SG_idhy); %So if the 1st Nov or after hydrological year is the same year
SG_T.HyYear(~SG_idhy)=SG_T.Year(~SG_idhy)-1; %So if before 1st Nov then hydrological year is the year before (note using HYear as year it starts in)

%Sums and averages
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_AD=grpstats(SG_T(:,2:end),'Datenum','mean'); %Calculates daily averages (NOTE THIS MIGHT NOT ALWAYS BE SUITABLE)
SG_SD=grpstats(SG_T(:,2:end),'Datenum','sum'); %Calculates daily sums
SG_AM=grpstats(SG_T(:,2:end),'Month','mean'); %Calculates monthly mean
SG_SM=grpstats(SG_T(:,2:end),'Month','sum'); %Calculates monthly sum
SG_AH_wet=grpstats(SG_T(SG_T.Season==1,2:end),'Hour','mean'); %Wet season hourly average (nice works!)
SG_AH_dry=grpstats(SG_T(SG_T.Season==0,2:end),'Hour','mean'); %Dry season hourly average 
SG_A_season=grpstats(SG_T(:,2:end),'Season','mean'); %Average for each season
SG_S_season=grpstats(SG_T(:,2:end),'Season','sum'); %Sum for each season
SG_AY=grpstats(SG_T(:,2:end),'Year','mean'); %Annual average
SG_SY=grpstats(SG_T(:,2:end),'Year','sum'); %Annual sum
SG_A_all=mean(SG_T{:,2:end}); %Mean of all columns
SG_S_all=sum(SG_T{:,2:end}); %Sum of all columns
SG_AHyY=grpstats(SG_T(:,2:end),'HyYear','mean'); %Hydrological year average
SG_SHyY=grpstats(SG_T(:,2:end),'HyYear','sum'); %Hydrological year sum

%For albedo
SG_minD=grpstats(SG_T(2:end,2:end),'Datenum','min'); %Calculates daily min, need 2: end as first value 0
SG_AMminD=grpstats(SG_minD(:,2:end),'min_Month','mean'); %So this takes the monthly average of the daily mean

%====== Artesonraju Glacier ================

AG_T=struct2table(AG); %Create table
%Do some extra calculations
AG_T.Tmelt=AG_T.Smelt+AG_T.Imelt; %Total melt 
AG_L= 1000.*(2501.3 - 2.361.*(AG_T.Ta)); %%% Latent heat vaporization/condensaition [J/kg]
AG_Ls= Lf+AG_L; %%% Latent Heat of sublimation  [J/kg]
AG_T.EICE_Wm2=(AG_T.EICE./3600).*AG_Ls;
AG_T.ESN_Wm2=(AG_T.ESN./3600).*AG_Ls;

%Sort time
AG_T.Year=year(AG_T.DateTime); %Calculate year
AG_T.Month=month(AG_T.DateTime); %Calculate month
AG_T.Day=day(AG_T.DateTime); %Calcualte day
AG_T.Hour=hour(AG_T.DateTime); %Calculate hour
AG_T.Datenum=datenum(dateshift(AG_T.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)
AG_Date=datetime(0001,AG_T.Month,AG_T.Day,AG_T.Hour,0,0);%This is the day and time but no year
AG_idhy=AG_Date>=HyYear_threshold; %So if after 1st Nov then = 1
AG_T.HyYear(AG_idhy)=AG_T.Year(AG_idhy); %So if the 1st Nov or after hydrological year is the same year
AG_T.HyYear(~AG_idhy)=AG_T.Year(~AG_idhy)-1; %So if before 1st Nov then hydrological year is the year before (note using HYear as year it starts in)

%Sums and averages
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_AD=grpstats(AG_T(:,2:end),'Datenum','mean'); %Calculates daily averages (NOTE THIS MIGHT NOT ALWAYS BE SUITABLE)
AG_SD=grpstats(AG_T(:,2:end),'Datenum','sum'); %Calculates daily sums
AG_AM=grpstats(AG_T(:,2:end),'Month','mean'); %Calculates monthly mean
AG_SM=grpstats(AG_T(:,2:end),'Month','sum'); %Calculates monthly sum
AG_AH_wet=grpstats(AG_T(AG_T.Season==1,2:end),'Hour','mean'); %Wet season hourly average (nice works!)
AG_AH_dry=grpstats(AG_T(AG_T.Season==0,2:end),'Hour','mean'); %Dry season hourly average 
AG_A_season=grpstats(AG_T(:,2:end),'Season','mean'); %Average for each season
AG_S_season=grpstats(AG_T(:,2:end),'Season','sum'); %Sum for each season
AG_AY=grpstats(AG_T(:,2:end),'Year','mean'); %Annual average
AG_SY=grpstats(AG_T(:,2:end),'Year','sum'); %Annual sum
AG_A_all=mean(AG_T{:,2:end}); %Mean of all columns
AG_S_all=sum(AG_T{:,2:end}); %Sum of all columns
AG_AHyY=grpstats(AG_T(:,2:end),'HyYear','mean'); %Hydrological year average
AG_SHyY=grpstats(AG_T(:,2:end),'HyYear','sum'); %Hydrological year sum

%For albedo
AG_minD=grpstats(AG_T(2:end,2:end),'Datenum','min'); %Calculates daily min, need 2: end as first value 0
AG_AMminD=grpstats(AG_minD(:,2:end),'min_Month','mean'); %So this takes the monthly average of the daily mean

%====== Cuchillacocha Glacier ================

CG_T=struct2table(CG); %Create table
%Do some extra calculations
CG_T.Tmelt=CG_T.Smelt+CG_T.Imelt; %Total melt
CG_L= 1000.*(2501.3 - 2.361.*(CG_T.Ta)); %%% Latent heat vaporization/condensaition [J/kg]
CG_Ls= Lf+CG_L; %%% Latent Heat of sublimation  [J/kg]
CG_T.EICE_Wm2=(CG_T.EICE./3600).*CG_Ls;
CG_T.ESN_Wm2=(CG_T.ESN./3600).*CG_Ls;

%Sort time
CG_T.Year=year(CG_T.DateTime); %Calculate year
CG_T.Month=month(CG_T.DateTime); %Calculate month
CG_T.Day=day(CG_T.DateTime); %Calcualte day
CG_T.Hour=hour(CG_T.DateTime); %Calculate hour
CG_T.Datenum=datenum(dateshift(CG_T.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)
CG_Date=datetime(0001,CG_T.Month,CG_T.Day,CG_T.Hour,0,0);%This is the day and time but no year
CG_idhy=CG_Date>=HyYear_threshold; %So if after 1st Nov then = 1
CG_T.HyYear(CG_idhy)=CG_T.Year(CG_idhy); %So if the 1st Nov or after hydrological year is the same year
CG_T.HyYear(~CG_idhy)=CG_T.Year(~CG_idhy)-1; %So if before 1st Nov then hydrological year is the year before (note using HYear as year it starts in)

%Sums and averages
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_AD=grpstats(CG_T(:,2:end),'Datenum','mean'); %Calculates daily averages (NOTE THIS MIGHT NOT ALWAYS BE SUITABLE)
CG_SD=grpstats(CG_T(:,2:end),'Datenum','sum'); %Calculates daily sums
CG_AM=grpstats(CG_T(:,2:end),'Month','mean'); %Calculates monthly mean
CG_SM=grpstats(CG_T(:,2:end),'Month','sum'); %Calculates monthly sum
CG_AH_wet=grpstats(CG_T(CG_T.Season==1,2:end),'Hour','mean'); %Wet season hourly average (nice works!)
CG_AH_dry=grpstats(CG_T(CG_T.Season==0,2:end),'Hour','mean'); %Dry season hourly average 
CG_A_season=grpstats(CG_T(:,2:end),'Season','mean'); %Average for each season
CG_S_season=grpstats(CG_T(:,2:end),'Season','sum'); %Sum for each season
CG_AY=grpstats(CG_T(:,2:end),'Year','mean'); %Annual average
CG_SY=grpstats(CG_T(:,2:end),'Year','sum'); %Annual sum
CG_A_all=mean(CG_T{:,2:end}); %Mean of all columns
CG_S_all=sum(CG_T{:,2:end}); %Sum of all columns
CG_AHyY=grpstats(CG_T(:,2:end),'HyYear','mean'); %Hydrological year average
CG_SHyY=grpstats(CG_T(:,2:end),'HyYear','sum'); %Hydrological year sum

CG_minD=grpstats(CG_T(2:end,2:end),'Datenum','min'); %Calculates daily min, need 2: end as first value 0
CG_AMminD=grpstats(CG_minD(:,2:end),'min_Month','mean'); %So this takes the monthly average of the daily mean

%====== Quelccaya Ice Cap ================

QIC_T=struct2table(QIC); %Create table
%Do some extra calculations
QIC_T.Tmelt=QIC_T.Smelt+QIC_T.Imelt; %Total melt
QIC_L= 1000.*(2501.3 - 2.361.*(QIC_T.Ta)); %%% Latent heat vaporization/condensaition [J/kg]
QIC_Ls= Lf+QIC_L; %%% Latent Heat of sublimation  [J/kg]
QIC_T.EICE_Wm2=(QIC_T.EICE./3600).*QIC_Ls;
QIC_T.ESN_Wm2=(QIC_T.ESN./3600).*QIC_Ls;

%Sort time
QIC_T.Year=year(QIC_T.DateTime); %Calculate year
QIC_T.Month=month(QIC_T.DateTime); %Calculate month
QIC_T.Day=day(QIC_T.DateTime); %Calcualte day
QIC_T.Hour=hour(QIC_T.DateTime); %Calculate hour
QIC_T.Datenum=datenum(dateshift(QIC_T.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)
QIC_Date=datetime(0001,QIC_T.Month,QIC_T.Day,QIC_T.Hour,0,0);%This is the day and time but no year
QIC_idhy=QIC_Date>=HyYear_threshold; %So if after 1st Nov then = 1
QIC_T.HyYear(QIC_idhy)=QIC_T.Year(QIC_idhy); %So if the 1st Nov or after hydrological year is the same year
QIC_T.HyYear(~QIC_idhy)=QIC_T.Year(~QIC_idhy)-1; %So if before 1st Nov then hydrological year is the year before (note using HYear as year it starts in)

%Sums and averages
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_AD=grpstats(QIC_T(:,2:end),'Datenum','mean'); %Calculates daily averages (NOTE THIS MIGHT NOT ALWAYS BE SUITABLE)
QIC_SD=grpstats(QIC_T(:,2:end),'Datenum','sum'); %Calculates daily sums
QIC_AM=grpstats(QIC_T(:,2:end),'Month','mean'); %Calculates monthly mean
QIC_SM=grpstats(QIC_T(:,2:end),'Month','sum'); %Calculates monthly sum
QIC_AH_wet=grpstats(QIC_T(QIC_T.Season==1,2:end),'Hour','mean'); %Wet season hourly average (nice works!)
QIC_AH_dry=grpstats(QIC_T(QIC_T.Season==0,2:end),'Hour','mean'); %Dry season hourly average 
QIC_A_season=grpstats(QIC_T(:,2:end),'Season','mean'); %Average for each season
QIC_S_season=grpstats(QIC_T(:,2:end),'Season','sum'); %Sum for each season
QIC_AY=grpstats(QIC_T(:,2:end),'Year','mean'); %Annual average
QIC_SY=grpstats(QIC_T(:,2:end),'Year','sum'); %Annual sum
QIC_A_all=mean(QIC_T{:,2:end}); %Mean of all columns
QIC_S_all=sum(QIC_T{:,2:end}); %Sum of all columns
QIC_AHyY=grpstats(QIC_T(:,2:end),'HyYear','mean'); %Hydrological year average
QIC_SHyY=grpstats(QIC_T(:,2:end),'HyYear','sum'); %Hydrological year sum

QIC_minD=grpstats(QIC_T(2:end,2:end),'Datenum','min'); %Calculates daily min, need 2: end as first value 0
QIC_AMminD=grpstats(QIC_minD(:,2:end),'min_Month','mean'); %So this takes the monthly average of the daily mean

%====== Quisoquipina Glacier ================

QQG_T=struct2table(QQG); %Create table
%Do some extra calculations
QQG_T.Tmelt=QQG_T.Smelt+QQG_T.Imelt; %Total melt
QQG_L= 1000.*(2501.3 - 2.361.*(QQG_T.Ta)); %%% Latent heat vaporization/condensaition [J/kg]
QQG_Ls= Lf+QQG_L; %%% Latent Heat of sublimation  [J/kg]
QQG_T.EICE_Wm2=(QQG_T.EICE./3600).*QQG_Ls;
QQG_T.ESN_Wm2=(QQG_T.ESN./3600).*QQG_Ls;

%Sort time
QQG_T.Year=year(QQG_T.DateTime); %Calculate year
QQG_T.Month=month(QQG_T.DateTime); %Calculate month
QQG_T.Day=day(QQG_T.DateTime); %Calcualte day
QQG_T.Hour=hour(QQG_T.DateTime); %Calculate hour
QQG_T.Datenum=datenum(dateshift(QQG_T.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)
QQG_Date=datetime(0001,QQG_T.Month,QQG_T.Day,QQG_T.Hour,0,0);%This is the day and time but no year
QQG_idhy=QQG_Date>=HyYear_threshold; %So if after 1st Nov then = 1
QQG_T.HyYear(QQG_idhy)=QQG_T.Year(QQG_idhy); %So if the 1st Nov or after hydrological year is the same year
QQG_T.HyYear(~QQG_idhy)=QQG_T.Year(~QQG_idhy)-1; %So if before 1st Nov then hydrological year is the year before (note using HYear as year it starts in)

%Sums and averages
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_AD=grpstats(QQG_T(:,2:end),'Datenum','mean'); %Calculates daily averages (NOTE THIS MIGHT NOT ALWAYS BE SUITABLE)
QQG_SD=grpstats(QQG_T(:,2:end),'Datenum','sum'); %Calculates daily sums
QQG_AM=grpstats(QQG_T(:,2:end),'Month','mean'); %Calculates monthly mean
QQG_SM=grpstats(QQG_T(:,2:end),'Month','sum'); %Calculates monthly sum
QQG_AH_wet=grpstats(QQG_T(QQG_T.Season==1,2:end),'Hour','mean'); %Wet season hourly average (nice works!)
QQG_AH_dry=grpstats(QQG_T(QQG_T.Season==0,2:end),'Hour','mean'); %Dry season hourly average 
QQG_A_season=grpstats(QQG_T(:,2:end),'Season','mean'); %Average for each season
QQG_S_season=grpstats(QQG_T(:,2:end),'Season','sum'); %Sum for each season
QQG_AY=grpstats(QQG_T(:,2:end),'Year','mean'); %Annual average
QQG_SY=grpstats(QQG_T(:,2:end),'Year','sum'); %Annual sum
QQG_A_all=mean(QQG_T{:,2:end}); %Mean of all columns
QQG_S_all=sum(QQG_T{:,2:end}); %Sum of all columns
QQG_AHyY=grpstats(QQG_T(:,2:end),'HyYear','mean'); %Hydrological year average
QQG_SHyY=grpstats(QQG_T(:,2:end),'HyYear','sum'); %Hydrological year sum

QQG_minD=grpstats(QQG_T(2:end,2:end),'Datenum','min'); %Calculates daily min, need 2: end as first value 0
QQG_AMminD=grpstats(QQG_minD(:,2:end),'min_Month','mean'); %So this takes the monthly average of the daily mean

%====== Right lots of graphs =================

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
%on;plot(Hours,All_SWnet); will plot with nice colours

Hours=[0:1:23]';

% 1. Big subplot of radiation AH
SG_AH_fluxes=SG_AH(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'});
SG_AH_wet_fluxes=SG_AH_wet(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'});
SG_AH_dry_fluxes=SG_AH_dry(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'});
AG_AH_fluxes=AG_AH(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'});
AG_AH_wet_fluxes=AG_AH_wet(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'});
AG_AH_dry_fluxes=AG_AH_dry(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'});
CG_AH_fluxes=CG_AH(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'}); 
CG_AH_wet_fluxes=CG_AH_wet(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'});
CG_AH_dry_fluxes=CG_AH_dry(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'});
QIC_AH_fluxes=QIC_AH(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'}); 
QIC_AH_wet_fluxes=QIC_AH_wet(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'});
QIC_AH_dry_fluxes=QIC_AH_dry(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'});
QQG_AH_fluxes=QQG_AH(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'});
QQG_AH_wet_fluxes=QQG_AH_wet(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'});
QQG_AH_dry_fluxes=QQG_AH_dry(:,{'mean_Rsw','mean_Rsw_out_mod','mean_Latm_mod','mean_L_out_mod','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G'});

%2. Big subplot of net radiation
SG_AH_netfluxes=SG_AH(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'});
SG_AH_wet_netfluxes=SG_AH_wet(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'});
SG_AH_dry_netfluxes=SG_AH_dry(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'});
AG_AH_netfluxes=AG_AH(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'});
AG_AH_wet_netfluxes=AG_AH_wet(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'});
AG_AH_dry_netfluxes=AG_AH_dry(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'});
CG_AH_netfluxes=CG_AH(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); 
CG_AH_wet_netfluxes=CG_AH_wet(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'});
CG_AH_dry_netfluxes=CG_AH_dry(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'});
QIC_AH_netfluxes=QIC_AH(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); 
QIC_AH_wet_netfluxes=QIC_AH_wet(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'});
QIC_AH_dry_netfluxes=QIC_AH_dry(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'});
QQG_AH_netfluxes=QQG_AH(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'});
QQG_AH_wet_netfluxes=QQG_AH_wet(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'});
QQG_AH_dry_netfluxes=QQG_AH_dry(:,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'});

%3. Comparing fluxes this time

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"]; %Note now swapped to QQG then QIC
All_SWnet=[SG_AH_netfluxes{:,1},AG_AH_netfluxes{:,1},CG_AH_netfluxes{:,1},QQG_AH_netfluxes{:,1},QIC_AH_netfluxes{:,1}];
All_LWnet=[SG_AH_netfluxes{:,2},AG_AH_netfluxes{:,2},CG_AH_netfluxes{:,2},QQG_AH_netfluxes{:,2},QIC_AH_netfluxes{:,2}];
All_H=[SG_AH_netfluxes{:,3},AG_AH_netfluxes{:,3},CG_AH_netfluxes{:,3},QQG_AH_netfluxes{:,3},QIC_AH_netfluxes{:,3}];
All_LE=[SG_AH_netfluxes{:,4},AG_AH_netfluxes{:,4},CG_AH_netfluxes{:,4},QQG_AH_netfluxes{:,4},QIC_AH_netfluxes{:,4}];
All_Qv=[SG_AH_netfluxes{:,5},AG_AH_netfluxes{:,5},CG_AH_netfluxes{:,5},QQG_AH_netfluxes{:,5},QIC_AH_netfluxes{:,5}];
All_Qfm=[SG_AH_netfluxes{:,6},AG_AH_netfluxes{:,6},CG_AH_netfluxes{:,6},QQG_AH_netfluxes{:,6},QIC_AH_netfluxes{:,6}];
All_G=[SG_AH_netfluxes{:,7},AG_AH_netfluxes{:,7},CG_AH_netfluxes{:,7},QQG_AH_netfluxes{:,7},QIC_AH_netfluxes{:,7}];
All_dQ=[SG_AH_netfluxes{:,8},AG_AH_netfluxes{:,8},CG_AH_netfluxes{:,8},QQG_AH_netfluxes{:,8},QIC_AH_netfluxes{:,8}];

All_SWnet_wet=[SG_AH_wet_netfluxes{:,1},AG_AH_wet_netfluxes{:,1},CG_AH_wet_netfluxes{:,1},QQG_AH_wet_netfluxes{:,1},QIC_AH_wet_netfluxes{:,1}];
All_LWnet_wet=[SG_AH_wet_netfluxes{:,2},AG_AH_wet_netfluxes{:,2},CG_AH_wet_netfluxes{:,2},QQG_AH_wet_netfluxes{:,2},QIC_AH_wet_netfluxes{:,2}];
All_H_wet=[SG_AH_wet_netfluxes{:,3},AG_AH_wet_netfluxes{:,3},CG_AH_wet_netfluxes{:,3},QQG_AH_wet_netfluxes{:,3},QIC_AH_wet_netfluxes{:,3}];
All_LE_wet=[SG_AH_wet_netfluxes{:,4},AG_AH_wet_netfluxes{:,4},CG_AH_wet_netfluxes{:,4},QQG_AH_wet_netfluxes{:,4},QIC_AH_wet_netfluxes{:,4}];
All_Qv_wet=[SG_AH_wet_netfluxes{:,5},AG_AH_wet_netfluxes{:,5},CG_AH_wet_netfluxes{:,5},QQG_AH_wet_netfluxes{:,5},QIC_AH_wet_netfluxes{:,5}];
All_Qfm_wet=[SG_AH_wet_netfluxes{:,6},AG_AH_wet_netfluxes{:,6},CG_AH_wet_netfluxes{:,6},QQG_AH_wet_netfluxes{:,6},QIC_AH_wet_netfluxes{:,6}];
All_G_wet=[SG_AH_wet_netfluxes{:,7},AG_AH_wet_netfluxes{:,7},CG_AH_wet_netfluxes{:,7},QQG_AH_wet_netfluxes{:,7},QIC_AH_wet_netfluxes{:,7}];
All_dQ_wet=[SG_AH_wet_netfluxes{:,8},AG_AH_wet_netfluxes{:,8},CG_AH_wet_netfluxes{:,8},QQG_AH_wet_netfluxes{:,8},QIC_AH_wet_netfluxes{:,8}];

All_SWnet_dry=[SG_AH_dry_netfluxes{:,1},AG_AH_dry_netfluxes{:,1},CG_AH_dry_netfluxes{:,1},QQG_AH_dry_netfluxes{:,1},QIC_AH_dry_netfluxes{:,1}];
All_LWnet_dry=[SG_AH_dry_netfluxes{:,2},AG_AH_dry_netfluxes{:,2},CG_AH_dry_netfluxes{:,2},QQG_AH_dry_netfluxes{:,2},QIC_AH_dry_netfluxes{:,2}];
All_H_dry=[SG_AH_dry_netfluxes{:,3},AG_AH_dry_netfluxes{:,3},CG_AH_dry_netfluxes{:,3},QQG_AH_dry_netfluxes{:,3},QIC_AH_dry_netfluxes{:,3}];
All_LE_dry=[SG_AH_dry_netfluxes{:,4},AG_AH_dry_netfluxes{:,4},CG_AH_dry_netfluxes{:,4},QQG_AH_dry_netfluxes{:,4},QIC_AH_dry_netfluxes{:,4}];
All_Qv_dry=[SG_AH_dry_netfluxes{:,5},AG_AH_dry_netfluxes{:,5},CG_AH_dry_netfluxes{:,5},QQG_AH_dry_netfluxes{:,5},QIC_AH_dry_netfluxes{:,5}];
All_Qfm_dry=[SG_AH_dry_netfluxes{:,6},AG_AH_dry_netfluxes{:,6},CG_AH_dry_netfluxes{:,6},QQG_AH_dry_netfluxes{:,6},QIC_AH_dry_netfluxes{:,6}];
All_G_dry=[SG_AH_dry_netfluxes{:,7},AG_AH_dry_netfluxes{:,7},CG_AH_dry_netfluxes{:,7},QQG_AH_dry_netfluxes{:,7},QIC_AH_dry_netfluxes{:,7}];
All_dQ_dry=[SG_AH_dry_netfluxes{:,8},AG_AH_dry_netfluxes{:,8},CG_AH_dry_netfluxes{:,8},QQG_AH_dry_netfluxes{:,8},QIC_AH_dry_netfluxes{:,8}];


All_SWnet_sdiff=All_SWnet_wet-All_SWnet_dry; %So wet - dry
All_LWnet_sdiff=All_LWnet_wet-All_LWnet_dry;
All_H_sdiff=All_H_wet-All_H_dry;
All_LE_sdiff=All_LE_wet-All_LE_dry;
All_Qv_sdiff=All_Qv_wet-All_Qv_dry;
All_Qfm_sdiff=All_Qfm_wet-All_Qfm_dry;
All_G_sdiff=All_G_wet-All_G_dry;
All_dQ_sdiff=All_dQ_wet-All_dQ_dry;

figure(3)
ax=subplot(8,4,1);ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_SWnet);
ylabel('S^\ast (W m^-^2)'); title ('S^\ast All');
xlim([0 23]);ylim([0 600]);
ax=subplot(8,4,2);ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_SWnet_wet);
title('S^\ast wet');
xlim([0 23]); ylim([0 600]);
ax=subplot(8,4,3); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_SWnet_dry);
title('S^\ast dry');
xlim([0 23]); ylim([0 600]);
ax=subplot(8,4,4);ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_SWnet_sdiff);
title('S^\ast wet-dry');
xlim([0 23]); %ylim([0 600]);


ax=subplot(8,4,5);ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_LWnet);
ylabel('L^\ast (W m^-^2)'); title ('L^\ast all');
xlim([0 23]); ylim([-70 20]);
ax=subplot(8,4,6); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_LWnet_wet);
title('L^\ast wet');
xlim([0 23]); ylim([-70 20]);
ax=subplot(8,4,7);ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_LWnet_dry);
title('L^\ast dry');
xlim([0 23]); ylim([-70 20]);
ax=subplot(8,4,8);ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_LWnet_sdiff);
title('L^\ast wet-dry');
xlim([0 23]); ylim([0 50]);

ax=subplot(8,4,9); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_H);
ylabel('H (W m^-^2)');title ('H All');
xlim([0 23]);ylim([-20 50]);
ax=subplot(8,4,10); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_H_wet);
title('H wet');
xlim([0 23]);ylim([-20 50]);
ax=subplot(8,4,11); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_H_dry);
title('H dry');
xlim([0 23]);ylim([-20 50]);
ax=subplot(8,4,12); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_H_sdiff);
title('H wet-dry');
xlim([0 23]);ylim([-20 20]);

ax=subplot(8,4,13);ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_LE);
ylabel('\lambdaE (W m^-^2)'); title ('\lambdaE All');
xlim([0 23]); ylim([-50 20]);
ax=subplot(8,4,14);ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_LE_wet);
title('\lambdaE wet');
xlim([0 23]); ylim([-50 20]);
ax=subplot(8,4,15);ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_LE_dry);
title('\lambdaE dry');
xlim([0 23]); ylim([-50 20]);
ax=subplot(8,4,16); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_LE_sdiff);
title('\lambdaE wet-dry');
xlim([0 23]); ylim([0 60]);

ax=subplot(8,4,17);ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Qv);
ylabel('Q_v (W m^-^2)');title ('Q_v All');
xlim([0 23]); ylim([-0.2 1]);
ax=subplot(8,4,18);ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Qv_wet);
title('Q_v wet');
xlim([0 23]);ylim([-0.2 1]);
ax=subplot(8,4,19);ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Qv_dry);
title('Q_v dry');
xlim([0 23]); ylim([-0.2 1]);
ax=subplot(8,4,20); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_Qv_sdiff);
title('Q_v wet-dry');
xlim([0 23]); ylim([-0.2 0.7]);

ax=subplot(8,4,21); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_Qfm);
ylabel('Q_f_m (W m^-^2)');title ('Q_f_m All');
xlim([0 23]); ylim([-6 6]);
ax=subplot(8,4,22);ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Qfm_wet);
title('Q_f_m wet');
xlim([0 23]);ylim([-6 6]);
ax=subplot(8,4,23); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_Qfm_dry);
title('Q_f_m dry');
xlim([0 23]); ylim([-6 6]);
ax=subplot(8,4,24); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_Qfm_sdiff);
title('Q_f_m wet-dry');
xlim([0 23]); %ylim([-6 6]);

ax=subplot(8,4,25); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_G);
ylabel('G (W m^-^2)');title ('G All');
xlim([0 23]); ylim([-2 4]);
ax=subplot(8,4,26); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_G_wet);
title('G wet');
xlim([0 23]);ylim([-2 4]);
ax=subplot(8,4,27); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_G_dry);
title('G dry');
xlim([0 23]); ylim([-2 4]);
ax=subplot(8,4,28); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_G_sdiff);
title('G wet-dry');
xlim([0 23]); %ylim([-2 4]);

ax=subplot(8,4,29); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_dQ);
ylabel('dQ (W m^-^2)');title ('dQ All');
xlim([0 23]); ylim([-100 550]);
ax=subplot(8,4,30); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_dQ_wet);
title('dQ wet');
xlim([0 23]);ylim([-100 550]);
ax=subplot(8,4,31); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_dQ_dry);
title('dQ dry');
xlim([0 23]); ylim([-100 550]);
ax=subplot(8,4,32); ax.ColorOrder = Fluxes_colours; hold on;plot(Hours,All_dQ_sdiff);
title('dQ wet-dry');
xlim([0 23]); %ylim([-2 4]);

f1L=legend(All_names,'Location','SouthOutside','Orientation','Horizontal');
f1P = [0.5 0.02 0.03 0.03];
set(f1L,'Position', f1P,'Units', 'normalized');

%4. Stacked bar charts on mean fluxes

SG_A_fluxes(1,:)=array2table(SG_A_all(:,[40 41 7 19 21 20 6 35])); %All
SG_A_fluxes(2,:)=SG_A_season(1,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); %Dry season
SG_A_fluxes(3,:)=SG_A_season(2,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); %Wet season
AG_A_fluxes(1,:)=array2table(AG_A_all(:,[40 41 7 19 21 20 6 35])); %All
AG_A_fluxes(2,:)=AG_A_season(1,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); %Dry season
AG_A_fluxes(3,:)=AG_A_season(2,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); %Wet season
CG_A_fluxes(1,:)=array2table(CG_A_all(:,[39 40 7 19 21 20 6 34])); %All %Columns in different order in CG
CG_A_fluxes(2,:)=CG_A_season(1,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); %Dry season
CG_A_fluxes(3,:)=CG_A_season(2,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); %Wet season
QIC_A_fluxes(1,:)=array2table(QIC_A_all(:,[39 40 7 19 21 20 6 34])); %All %Columns in different order in QIC
QIC_A_fluxes(2,:)=QIC_A_season(1,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); %Dry season
QIC_A_fluxes(3,:)=QIC_A_season(2,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); %Wet season
QQG_A_fluxes(1,:)=array2table(QQG_A_all(:,[40 41 7 19 21 20 6 35])); %All
QQG_A_fluxes(2,:)=QQG_A_season(1,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); %Dry season
QQG_A_fluxes(3,:)=QQG_A_season(2,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); %Wet season

Gap=zeros(1,8);
All_A_fluxes=[table2array(SG_A_fluxes);Gap;table2array(AG_A_fluxes);Gap;table2array(CG_A_fluxes);Gap;table2array(QQG_A_fluxes);Gap;table2array(QIC_A_fluxes)];
All_A_fluxes(4,8)=NaN; All_A_fluxes(8,8)=NaN; All_A_fluxes(12,8)=NaN; All_A_fluxes(16,8)=NaN; %Add NaNs in dQ

%figure(4)
SG_labels=categorical({'SG All','SG Dry','SG Wet'}); %This is the notation for categorical data
% bar(SG_labels,SG_A_fluxes,'stacked');
% legend('SW net','LW net','H','LE','Qv','Qfm','G');

figure5=figure(5);
All_labels=categorical({'All','Dry','Wet','All','Dry','Wet','All','Dry','Wet','All','Dry','Wet','All','Dry','Wet'});
ax=axes; ax.ColorOrder = Fluxes_colours; hold on; bar(All_A_fluxes(:,1:7),'stacked'); hold on; %Not including dQ
scatter([1:19],All_A_fluxes(:,8),22,'w','d','filled','MarkerEdgeColor',Fluxes_colours(1,:));
legend('S^\ast','L^\ast','H','\lambdaE','Q_v','Q_f_m','G','dQ','Location','SouthOutside','Orientation','Horizontal');
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('Energy fluxes (W m^-^2)'); 
% Create textbox
annotation(figure5,'textbox',...
    [0.181293956776012 0.861818658280922 0.0638967551622419 0.0534591194968551],'String','SG','HorizontalAlignment','center','FitBoxToText','off','EdgeColor',[1 1 1]);
annotation(figure5,'textbox',...
    [0.332939438421634 0.859329140461216 0.0638967551622417 0.0534591194968551],'String','AG','HorizontalAlignment','center','FitBoxToText','off','EdgeColor',[1 1 1]);
annotation(figure5,'textbox',...
    [0.485377313094198 0.859800838574423 0.0638967551622417 0.0534591194968551],'String','CG','HorizontalAlignment','center','FitBoxToText','off','EdgeColor',[1 1 1]);
annotation(figure5,'textbox',...
    [0.637184857423796 0.860167714884695 0.0638967551622419 0.0534591194968551],'String','QQG','HorizontalAlignment','center','FitBoxToText','off','EdgeColor',[1 1 1]);
annotation(figure5,'textbox',...
    [0.788992401753394 0.858228511530398 0.0638967551622417 0.0534591194968551],'String','QIC','HorizontalAlignment','center','FitBoxToText','off','EdgeColor',[1 1 1]);

%6. Stacked bar chart on monthly mean fluxes for all stations
SG_AM_fluxes=[SG_AM.mean_SWnet SG_AM.mean_LWnet SG_AM.mean_H SG_AM.mean_QE SG_AM.mean_Qv SG_AM.mean_Qfm SG_AM.mean_G SG_AM.mean_dQ];
AG_AM_fluxes=[AG_AM.mean_SWnet AG_AM.mean_LWnet AG_AM.mean_H AG_AM.mean_QE AG_AM.mean_Qv AG_AM.mean_Qfm AG_AM.mean_G AG_AM.mean_dQ];
CG_AM_fluxes=[CG_AM.mean_SWnet CG_AM.mean_LWnet CG_AM.mean_H CG_AM.mean_QE CG_AM.mean_Qv CG_AM.mean_Qfm CG_AM.mean_G CG_AM.mean_dQ];
QIC_AM_fluxes=[QIC_AM.mean_SWnet QIC_AM.mean_LWnet QIC_AM.mean_H QIC_AM.mean_QE QIC_AM.mean_Qv QIC_AM.mean_Qfm QIC_AM.mean_G QIC_AM.mean_dQ];
QQG_AM_fluxes=[QQG_AM.mean_SWnet QQG_AM.mean_LWnet QQG_AM.mean_H QQG_AM.mean_QE QQG_AM.mean_Qv QQG_AM.mean_Qfm QQG_AM.mean_G QQG_AM.mean_dQ];


figure(6)
Month_labels=categorical({'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'});
Month_labels = reordercats(Month_labels,{'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'});

ax=subplot(5,1,1); ax.ColorOrder = Fluxes_colours; hold on; bar(Month_labels,SG_AM_fluxes(:,1:7),'stacked'); hold on;
scatter(Month_labels,SG_AM_fluxes(:,8),22,'w','d','filled','MarkerEdgeColor',Fluxes_colours(1,:));
title('Shallap Glacier (4790 m)'); ylim([-120 200]);

ax=subplot(5,1,2);ax.ColorOrder = Fluxes_colours; hold on; bar(Month_labels,AG_AM_fluxes(:,1:7),'stacked');
scatter(Month_labels,AG_AM_fluxes(:,8),22,'w','d','filled','MarkerEdgeColor',Fluxes_colours(1,:));
title('Artesonraju Glacier (4797 m)');ylim([-120 200]);

ax=subplot(5,1,3);ax.ColorOrder = Fluxes_colours; hold on; bar(Month_labels,CG_AM_fluxes(:,1:7),'stacked');
scatter(Month_labels,CG_AM_fluxes(:,8),22,'w','d','filled','MarkerEdgeColor',Fluxes_colours(1,:));
title('Cuchillacocha Glacier (4821 m)'); ylabel('Energy fluxes (W m^-^2)');ylim([-120 200]);

ax=subplot(5,1,5);ax.ColorOrder = Fluxes_colours; hold on; bar(Month_labels,QIC_AM_fluxes(:,1:7),'stacked');
scatter(Month_labels,QIC_AM_fluxes(:,8),22,'w','d','filled','MarkerEdgeColor',Fluxes_colours(1,:));
title('Quelccaya Ice Cap (5650 m)'); ylim([-120 200]);

ax=subplot(5,1,4);ax.ColorOrder = Fluxes_colours; hold on; bar(Month_labels,QQG_AM_fluxes(:,1:7),'stacked');
scatter(Month_labels,QQG_AM_fluxes(:,8),22,'w','d','filled','MarkerEdgeColor',Fluxes_colours(1,:));
title('Quisoquipina Glacier (5180 m)');ylim([-120 200]);

Flux_names=["S^\ast","L^\ast","H","\lambdaE","Q_v","Q_f_m","G","dQ"];
f1L=legend(Flux_names,'Location','SouthOutside','Orientation','Horizontal');
f1P = [0.5 0.02 0.03 0.03];
set(f1L,'Position', f1P,'Units', 'normalized');

% ==== Weather conditions, so Ta, WS, RH, Swin and Lwin, Pr, Alb?  =====

%7. AH weather conditions
SG_AH_meteo=[SG_AH.mean_Ta SG_AH.mean_Ws SG_AH.mean_U SG_AH.mean_Rsw SG_AH.mean_Latm_mod SG_AH.mean_ALB SG_AH.mean_Pr];
AG_AH_meteo=[AG_AH.mean_Ta AG_AH.mean_Ws AG_AH.mean_U AG_AH.mean_Rsw AG_AH.mean_Latm_mod AG_AH.mean_ALB AG_AH.mean_Pr];
CG_AH_meteo=[CG_AH.mean_Ta CG_AH.mean_Ws CG_AH.mean_U CG_AH.mean_Rsw CG_AH.mean_Latm_mod CG_AH.mean_ALB CG_AH.mean_Pr];
QIC_AH_meteo=[QIC_AH.mean_Ta QIC_AH.mean_Ws QIC_AH.mean_U QIC_AH.mean_Rsw QIC_AH.mean_Latm_mod QIC_AH.mean_ALB QIC_AH.mean_Pr];
QQG_AH_meteo=[QQG_AH.mean_Ta QQG_AH.mean_Ws QQG_AH.mean_U QQG_AH.mean_Rsw QQG_AH.mean_Latm_mod QQG_AH.mean_ALB QQG_AH.mean_Pr];

SG_AH_wetmeteo=[SG_AH_wet.mean_Ta SG_AH_wet.mean_Ws SG_AH_wet.mean_U SG_AH_wet.mean_Rsw SG_AH_wet.mean_Latm_mod SG_AH_wet.mean_ALB SG_AH_wet.mean_Pr];
AG_AH_wetmeteo=[AG_AH_wet.mean_Ta AG_AH_wet.mean_Ws AG_AH_wet.mean_U AG_AH_wet.mean_Rsw AG_AH_wet.mean_Latm_mod AG_AH_wet.mean_ALB AG_AH_wet.mean_Pr];
CG_AH_wetmeteo=[CG_AH_wet.mean_Ta CG_AH_wet.mean_Ws CG_AH_wet.mean_U CG_AH_wet.mean_Rsw CG_AH_wet.mean_Latm_mod CG_AH_wet.mean_ALB CG_AH_wet.mean_Pr];
QIC_AH_wetmeteo=[QIC_AH_wet.mean_Ta QIC_AH_wet.mean_Ws QIC_AH_wet.mean_U QIC_AH_wet.mean_Rsw QIC_AH_wet.mean_Latm_mod QIC_AH_wet.mean_ALB QIC_AH_wet.mean_Pr];
QQG_AH_wetmeteo=[QQG_AH_wet.mean_Ta QQG_AH_wet.mean_Ws QQG_AH_wet.mean_U QQG_AH_wet.mean_Rsw QQG_AH_wet.mean_Latm_mod QQG_AH_wet.mean_ALB QQG_AH_wet.mean_Pr];

SG_AH_drymeteo=[SG_AH_dry.mean_Ta SG_AH_dry.mean_Ws SG_AH_dry.mean_U SG_AH_dry.mean_Rsw SG_AH_dry.mean_Latm_mod SG_AH_dry.mean_ALB SG_AH_dry.mean_Pr];
AG_AH_drymeteo=[AG_AH_dry.mean_Ta AG_AH_dry.mean_Ws AG_AH_dry.mean_U AG_AH_dry.mean_Rsw AG_AH_dry.mean_Latm_mod AG_AH_dry.mean_ALB AG_AH_dry.mean_Pr];
CG_AH_drymeteo=[CG_AH_dry.mean_Ta CG_AH_dry.mean_Ws CG_AH_dry.mean_U CG_AH_dry.mean_Rsw CG_AH_dry.mean_Latm_mod CG_AH_dry.mean_ALB CG_AH_dry.mean_Pr];
QIC_AH_drymeteo=[QIC_AH_dry.mean_Ta QIC_AH_dry.mean_Ws QIC_AH_dry.mean_U QIC_AH_dry.mean_Rsw QIC_AH_dry.mean_Latm_mod QIC_AH_dry.mean_ALB QIC_AH_dry.mean_Pr];
QQG_AH_drymeteo=[QQG_AH_dry.mean_Ta QQG_AH_dry.mean_Ws QQG_AH_dry.mean_U QQG_AH_dry.mean_Rsw QQG_AH_dry.mean_Latm_mod QQG_AH_dry.mean_ALB QQG_AH_dry.mean_Pr];

%I've swapped QQG and QIC round now
All_Ta_AH=[SG_AH_meteo(:,1),AG_AH_meteo(:,1),CG_AH_meteo(:,1),QQG_AH_meteo(:,1),QIC_AH_meteo(:,1)];
All_Ws_AH=[SG_AH_meteo(:,2),AG_AH_meteo(:,2),CG_AH_meteo(:,2),QQG_AH_meteo(:,2),QIC_AH_meteo(:,2)];
All_U_AH=[SG_AH_meteo(:,3),AG_AH_meteo(:,3),CG_AH_meteo(:,3),QQG_AH_meteo(:,3),QIC_AH_meteo(:,3)];
All_Rsw_AH=[SG_AH_meteo(:,4),AG_AH_meteo(:,4),CG_AH_meteo(:,4),QQG_AH_meteo(:,4),QIC_AH_meteo(:,4)];
All_Latm_mod_AH=[SG_AH_meteo(:,5),AG_AH_meteo(:,5),CG_AH_meteo(:,5),QQG_AH_meteo(:,5),QIC_AH_meteo(:,5)];
All_ALB_AH=[SG_AH_meteo(:,6),AG_AH_meteo(:,6),CG_AH_meteo(:,6),QQG_AH_meteo(:,6),QIC_AH_meteo(:,6)];
All_Pr_AH=[SG_AH_meteo(:,7),AG_AH_meteo(:,7),CG_AH_meteo(:,7),QQG_AH_meteo(:,7),QIC_AH_meteo(:,7)];

All_Ta_wetAH=[SG_AH_wetmeteo(:,1),AG_AH_wetmeteo(:,1),CG_AH_wetmeteo(:,1),QQG_AH_wetmeteo(:,1),QIC_AH_wetmeteo(:,1)];
All_Ws_wetAH=[SG_AH_wetmeteo(:,2),AG_AH_wetmeteo(:,2),CG_AH_wetmeteo(:,2),QQG_AH_wetmeteo(:,2),QIC_AH_wetmeteo(:,2)];
All_U_wetAH=[SG_AH_wetmeteo(:,3),AG_AH_wetmeteo(:,3),CG_AH_wetmeteo(:,3),QQG_AH_wetmeteo(:,3),QIC_AH_wetmeteo(:,3)];
All_Rsw_wetAH=[SG_AH_wetmeteo(:,4),AG_AH_wetmeteo(:,4),CG_AH_wetmeteo(:,4),QQG_AH_wetmeteo(:,4),QIC_AH_wetmeteo(:,4)];
All_Latm_mod_wetAH=[SG_AH_wetmeteo(:,5),AG_AH_wetmeteo(:,5),CG_AH_wetmeteo(:,5),QQG_AH_wetmeteo(:,5),QIC_AH_wetmeteo(:,5)];
All_ALB_wetAH=[SG_AH_wetmeteo(:,6),AG_AH_wetmeteo(:,6),CG_AH_wetmeteo(:,6),QQG_AH_wetmeteo(:,6),QIC_AH_wetmeteo(:,6)];
All_Pr_wetAH=[SG_AH_wetmeteo(:,7),AG_AH_wetmeteo(:,7),CG_AH_wetmeteo(:,7),QQG_AH_wetmeteo(:,7),QIC_AH_wetmeteo(:,7)];

All_Ta_dryAH=[SG_AH_drymeteo(:,1),AG_AH_drymeteo(:,1),CG_AH_drymeteo(:,1),QQG_AH_drymeteo(:,1),QIC_AH_drymeteo(:,1)];
All_Ws_dryAH=[SG_AH_drymeteo(:,2),AG_AH_drymeteo(:,2),CG_AH_drymeteo(:,2),QQG_AH_drymeteo(:,2),QIC_AH_drymeteo(:,2)];
All_U_dryAH=[SG_AH_drymeteo(:,3),AG_AH_drymeteo(:,3),CG_AH_drymeteo(:,3),QQG_AH_drymeteo(:,3),QIC_AH_drymeteo(:,3)];
All_Rsw_dryAH=[SG_AH_drymeteo(:,4),AG_AH_drymeteo(:,4),CG_AH_drymeteo(:,4),QQG_AH_drymeteo(:,4),QIC_AH_drymeteo(:,4)];
All_Latm_mod_dryAH=[SG_AH_drymeteo(:,5),AG_AH_drymeteo(:,5),CG_AH_drymeteo(:,5),QQG_AH_drymeteo(:,5),QIC_AH_drymeteo(:,5)];
All_ALB_dryAH=[SG_AH_drymeteo(:,6),AG_AH_drymeteo(:,6),CG_AH_drymeteo(:,6),QQG_AH_drymeteo(:,6),QIC_AH_drymeteo(:,6)];
All_Pr_dryAH=[SG_AH_drymeteo(:,7),AG_AH_drymeteo(:,7),CG_AH_drymeteo(:,7),QQG_AH_drymeteo(:,7),QIC_AH_drymeteo(:,7)];

All_Ta_AH_sdiff=All_Ta_wetAH-All_Ta_dryAH;
All_Ws_AH_sdiff=All_Ws_wetAH-All_Ws_dryAH;
All_U_AH_sdiff=All_U_wetAH-All_U_dryAH;
All_Rsw_AH_sdiff=All_Rsw_wetAH-All_Rsw_dryAH;
All_Latm_mod_AH_sdiff=All_Latm_mod_wetAH-All_Latm_mod_dryAH;
All_ALB_AH_sdiff=All_ALB_wetAH-All_ALB_dryAH;
All_Pr_AH_sdiff=All_Pr_wetAH-All_Pr_dryAH;


figure(7)
ax=subplot(6,4,1); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Ta_AH);
ylabel('T_a ({\circ}C)'); title ('T_a All');
xlim([0 23]);ylim([-10 5]);
ax=subplot(6,4,2); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Ta_wetAH); title('T_a wet');
xlim([0 23]); ylim([-10 5]);
ax=subplot(6,4,3); ax.ColorOrder = Fluxes_colours;hold on;  plot(Hours,All_Ta_dryAH); title('T_a dry');
xlim([0 23]); ylim([-10 5]);
ax=subplot(6,4,4); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Ta_AH_sdiff); title('T_a wet-dry');
xlim([0 23]); ylim([-1.5 3]);

ax=subplot(6,4,5);ax.ColorOrder = Fluxes_colours;hold on; plot(Hours,All_Ws_AH);
ylabel('WS (m s^-^1)'); title ('WS all');
xlim([0 23]); ylim([0 8]);
ax=subplot(6,4,6); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Ws_wetAH);
title('WS wet');
xlim([0 23]); ylim([0 8]);
ax=subplot(6,4,7); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Ws_dryAH);
title('WS dry');
xlim([0 23]); ylim([0 8]);
ax=subplot(6,4,8); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Ws_AH_sdiff);
title('WS wet-dry');
xlim([0 23]); ylim([-2.5 1]);

ax=subplot(6,4,9); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_U_AH);
ylabel('RH (fraction)');title ('RH All');
xlim([0 23]);ylim([0.5 1]);
ax=subplot(6,4,10); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_U_wetAH);
title('RH wet');
xlim([0 23]);ylim([0.5 1]);
ax=subplot(6,4,11); ax.ColorOrder = Fluxes_colours;hold on;  plot(Hours,All_U_dryAH);
title('RH dry');
xlim([0 23]);ylim([0.5 1]);
ax=subplot(6,4,12); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_U_AH_sdiff);
title('RH wet-dry');
xlim([0 23]);ylim([0.05 0.25]);

ax=subplot(6,4,13);ax.ColorOrder = Fluxes_colours;hold on;  plot(Hours,All_Rsw_AH);
ylabel('S\downarrow (W m^-^2)'); title ('S\downarrow All');
xlim([0 23]); ylim([0 1000]);
ax=subplot(6,4,14);ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Rsw_wetAH);
title('S\downarrow wet');
xlim([0 23]); ylim([0 1000]);
ax=subplot(6,4,15); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Rsw_dryAH);
title('S\downarrow dry');
xlim([0 23]); ylim([0 1000]);
ax=subplot(6,4,16); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Rsw_AH_sdiff);
title('S\downarrow wet-dry');
xlim([0 23]); %ylim([0 1000]);

ax=subplot(6,4,17); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Latm_mod_AH);
ylabel('L\downarrow (W m^-^2)');title ('L\downarrow All');
xlim([0 23]); ylim([200 330]);
ax=subplot(6,4,18); ax.ColorOrder = Fluxes_colours;hold on;  plot(Hours,All_Latm_mod_wetAH);
title('L\downarrow wet');
xlim([0 23]);ylim([200 330]);
ax=subplot(6,4,19); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Latm_mod_dryAH);
title('L\downarrow dry');
xlim([0 23]); ylim([200 330]);
ax=subplot(6,4,20); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Latm_mod_AH_sdiff);
title('L\downarrow wet-dry');
xlim([0 23]); ylim([10 60]);

ax=subplot(6,4,21); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Pr_AH);
ylabel('P_r (mm h^-^1)');title ('P_r All');
xlim([0 23]); ylim([0 1.3]);
ax=subplot(6,4,22); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Pr_wetAH);
title('P_r wet');
xlim([0 23]);ylim([0 1.3]);
ax=subplot(6,4,23); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Pr_dryAH);
title('P_r dry');
xlim([0 23]); ylim([0 1.3]);
ax=subplot(6,4,24); ax.ColorOrder = Fluxes_colours; hold on; plot(Hours,All_Pr_AH_sdiff);
title('P_r wet-dry');
xlim([0 23]); ylim([0 1.3]);

f1L=legend(All_names,'Location','SouthOutside','Orientation','Horizontal');
f1P = [0.5 0.02 0.03 0.03];
set(f1L,'Position', f1P,'Units', 'normalized');

%Monthly weather conditions
%I've swapped QQG and QIC
All_Ta_AM=[SG_AM.mean_Ta,AG_AM.mean_Ta,CG_AM.mean_Ta,QQG_AM.mean_Ta,QIC_AM.mean_Ta];
All_WS_AM=[SG_AM.mean_Ws,AG_AM.mean_Ws,CG_AM.mean_Ws,QQG_AM.mean_Ws,QIC_AM.mean_Ws];
All_U_AM=[SG_AM.mean_U,AG_AM.mean_U,CG_AM.mean_U,QQG_AM.mean_U,QIC_AM.mean_U];
All_Rsw_AM=[SG_AM.mean_Rsw,AG_AM.mean_Rsw,CG_AM.mean_Rsw,QQG_AM.mean_Rsw,QIC_AM.mean_Rsw];
All_Latm_mod_AM=[SG_AM.mean_Latm_mod,AG_AM.mean_Latm_mod,CG_AM.mean_Latm_mod,QQG_AM.mean_Latm_mod,QIC_AM.mean_Latm_mod];
All_ALB_AM=[SG_AM.mean_ALB,AG_AM.mean_ALB,CG_AM.mean_ALB,QQG_AM.mean_ALB,QIC_AM.mean_ALB];
All_Pr_SM=[SG_SM.sum_Pr,AG_SM.sum_Pr,CG_SM.sum_Pr,QQG_SM.sum_Pr,QIC_SM.sum_Pr]; %This is the sum of precipitation

%========= Snow and ice melt and snow depth ================

SG_SM_melt=[SG_SM.sum_Imelt SG_SM.sum_Smelt SG_SM.sum_Tmelt]; %NOTE THESE ARE SUMS OF ALL THE TIME!! - so not average monthly sum but just the actual sum
AG_SM_melt=[AG_SM.sum_Imelt AG_SM.sum_Smelt AG_SM.sum_Tmelt];
CG_SM_melt=[CG_SM.sum_Imelt CG_SM.sum_Smelt CG_SM.sum_Tmelt];
QIC_SM_melt=[QIC_SM.sum_Imelt QIC_SM.sum_Smelt QIC_SM.sum_Tmelt];
QQG_SM_melt=[QQG_SM.sum_Imelt QQG_SM.sum_Smelt QQG_SM.sum_Tmelt];

SG_AM_melt=[SG_AM.mean_Imelt SG_AM.mean_Smelt SG_AM.mean_Tmelt]; %Mean melt this time
AG_AM_melt=[AG_AM.mean_Imelt AG_AM.mean_Smelt AG_AM.mean_Tmelt];
CG_AM_melt=[CG_AM.mean_Imelt CG_AM.mean_Smelt CG_AM.mean_Tmelt];
QIC_AM_melt=[QIC_AM.mean_Imelt QIC_AM.mean_Smelt QIC_AM.mean_Tmelt];
QQG_AM_melt=[QQG_AM.mean_Imelt QQG_AM.mean_Smelt QQG_AM.mean_Tmelt];

SG_AM_abla=[SG_AM.mean_Imelt SG_AM.mean_Smelt SG_AM.mean_Tmelt SG_AM.mean_EICE SG_AM.mean_ESN SG_AM.mean_EICE+SG_AM.mean_ESN]; %Mean ablation, so including sublimation, last column is sum of sublimation from ice and snow
AG_AM_abla=[AG_AM.mean_Imelt AG_AM.mean_Smelt AG_AM.mean_Tmelt AG_AM.mean_EICE AG_AM.mean_ESN AG_AM.mean_EICE+AG_AM.mean_ESN];
CG_AM_abla=[CG_AM.mean_Imelt CG_AM.mean_Smelt CG_AM.mean_Tmelt CG_AM.mean_EICE CG_AM.mean_ESN CG_AM.mean_EICE+CG_AM.mean_ESN];
QIC_AM_abla=[QIC_AM.mean_Imelt QIC_AM.mean_Smelt QIC_AM.mean_Tmelt QIC_AM.mean_EICE QIC_AM.mean_ESN QIC_AM.mean_EICE+QIC_AM.mean_ESN];
QQG_AM_abla=[QQG_AM.mean_Imelt QQG_AM.mean_Smelt QQG_AM.mean_Tmelt QQG_AM.mean_EICE QQG_AM.mean_ESN QQG_AM.mean_EICE+QQG_AM.mean_ESN];

%Perhaps also include %snow melt of total ablation

SG_SM.PercSmelt=(SG_SM.sum_Smelt./(SG_SM.sum_EICE+SG_SM.sum_ESN+SG_SM.sum_Imelt+SG_SM.sum_Smelt)).*100;
AG_SM.PercSmelt=(AG_SM.sum_Smelt./(AG_SM.sum_EICE+AG_SM.sum_ESN+AG_SM.sum_Imelt+AG_SM.sum_Smelt)).*100;
CG_SM.PercSmelt=(CG_SM.sum_Smelt./(CG_SM.sum_EICE+CG_SM.sum_ESN+CG_SM.sum_Imelt+CG_SM.sum_Smelt)).*100;
QIC_SM.PercSmelt=(QIC_SM.sum_Smelt./(QIC_SM.sum_EICE+QIC_SM.sum_ESN+QIC_SM.sum_Imelt+QIC_SM.sum_Smelt)).*100;
QQG_SM.PercSmelt=(QQG_SM.sum_Smelt./(QQG_SM.sum_EICE+QQG_SM.sum_ESN+QQG_SM.sum_Imelt+QQG_SM.sum_Smelt)).*100;

%And perhaps % sublimation of total ablation
SG_SM.PercSub=((SG_SM.sum_EICE+SG_SM.sum_ESN)./(SG_SM.sum_EICE+SG_SM.sum_ESN+SG_SM.sum_Imelt+SG_SM.sum_Smelt)).*100;
AG_SM.PercSub=((AG_SM.sum_EICE+AG_SM.sum_ESN)./(AG_SM.sum_EICE+AG_SM.sum_ESN+AG_SM.sum_Imelt+AG_SM.sum_Smelt)).*100;
CG_SM.PercSub=((CG_SM.sum_EICE+CG_SM.sum_ESN)./(CG_SM.sum_EICE+CG_SM.sum_ESN+CG_SM.sum_Imelt+CG_SM.sum_Smelt)).*100;
QIC_SM.PercSub=((QIC_SM.sum_EICE+QIC_SM.sum_ESN)./(QIC_SM.sum_EICE+QIC_SM.sum_ESN+QIC_SM.sum_Imelt+QIC_SM.sum_Smelt)).*100;
QQG_SM.PercSub=((QQG_SM.sum_EICE+QQG_SM.sum_ESN)./(QQG_SM.sum_EICE+QQG_SM.sum_ESN+QQG_SM.sum_Imelt+QQG_SM.sum_Smelt)).*100;

%Melt, SND, % snow melt and sublimation and albedo
figure();
subplot(5,1,1); ba = bar(Month_labels,SG_AM_abla(:,[1:2 6]),'FaceColor','flat'); hold on; plot(Month_labels,SG_AMminD.mean_min_ALB,'Color','k','LineStyle','-.','LineWidth',0.7);
plot(Month_labels,SG_SM.PercSmelt./100,'Color',[0.6 0.6 0.6],'LineStyle',':','LineWidth',1); 
plot(Month_labels,SG_SM.PercSub./100,'Color',[0.9 0.17 0.71],'LineStyle',':','LineWidth',1);
ba(1).CData = [0.49 0.73 0.91]; ba(2).CData = [1 1 1]*0.8; ba(3).CData = [0.95 0.61 0.73];
title('Shallap Glacier (4790 m)'); ylim([0 1]);
yyaxis right; plot(Month_labels,SG_AM.mean_SND,'Color',[0.05 0.6 0.73],'LineStyle','--','LineWidth',1); hold on; set(gca,'YColor',[0.05 0.6 0.73]); ylim([0 0.7]);

subplot(5,1,2);ba2 = bar(Month_labels,AG_AM_abla(:,[1:2 6]),'FaceColor','flat');hold on; plot(Month_labels,AG_AMminD.mean_min_ALB,'Color','k','LineStyle','-.','LineWidth',0.7);
plot(Month_labels,AG_SM.PercSmelt./100,'Color',[0.6 0.6 0.6],'LineStyle',':','LineWidth',1);
plot(Month_labels,AG_SM.PercSub./100,'Color',[0.9 0.17 0.71],'LineStyle',':','LineWidth',1);
ba2(1).CData = [0.49 0.73 0.91]; ba2(2).CData = [1 1 1]*0.8; ba2(3).CData = [0.95 0.61 0.73];
title('Artesonraju Glacier (4797 m)');ylim([0 1]);
yyaxis right; plot(Month_labels,AG_AM.mean_SND,'Color',[0.05 0.6 0.73],'LineStyle','--','LineWidth',1); set(gca,'YColor',[0.05 0.6 0.73]);ylim([0 0.7]);

subplot(5,1,3);ba3 = bar(Month_labels,CG_AM_abla(:,[1:2 6]),'FaceColor','flat');hold on; plot(Month_labels,CG_AMminD.mean_min_ALB,'Color','k','LineStyle','-.','LineWidth',0.7);
plot(Month_labels,CG_SM.PercSmelt./100,'Color',[0.6 0.6 0.6],'LineStyle',':','LineWidth',1);
plot(Month_labels,CG_SM.PercSub./100,'Color',[0.9 0.17 0.71],'LineStyle',':','LineWidth',1);
ba3(1).CData = [0.49 0.73 0.91]; ba3(2).CData = [1 1 1]*0.8; ba3(3).CData = [0.95 0.61 0.73];
title('Cuchillacocha Glacier (4821 m)'); ylabel('Mean ablation (mm w.e. h^-^1) / Fraction snowmelt / Fraction sublimation / Albedo','FontSize',12);ylim([0 1]);
yyaxis right; plot(Month_labels,CG_AM.mean_SND,'Color',[0.05 0.6 0.73],'LineStyle','--','LineWidth',1); ylabel('Mean snow depth (m)','Color',[0.05 0.6 0.73],'FontSize',12); set(gca,'YColor',[0.05 0.6 0.73]);ylim([0 0.7]);

subplot(5,1,5);ba4 = bar(Month_labels,QIC_AM_abla(:,[1:2 6]),'FaceColor','flat');hold on; plot(Month_labels,QIC_AMminD.mean_min_ALB,'Color','k','LineStyle','-.','LineWidth',0.7);
plot(Month_labels,QIC_SM.PercSmelt./100,'Color',[0.6 0.6 0.6],'LineStyle',':','LineWidth',1);
plot(Month_labels,QIC_SM.PercSub./100,'Color',[0.9 0.17 0.71],'LineStyle',':','LineWidth',1);
ba4(1).CData = [0.49 0.73 0.91]; ba4(2).CData = [1 1 1]*0.8; ba4(3).CData = [0.95 0.61 0.73];
title('Quelccaya Ice Cap (5650 m)');ylim([0 1]);
yyaxis right; plot(Month_labels,QIC_AM.mean_SND,'Color',[0.05 0.6 0.73],'LineStyle','--','LineWidth',1); set(gca,'YColor',[0.05 0.6 0.73]);

subplot(5,1,4);ba5 = bar(Month_labels,QQG_AM_abla(:,[1:2 6]),'FaceColor','flat');hold on; plot(Month_labels,QQG_AMminD.mean_min_ALB,'Color','k','LineStyle','-.','LineWidth',0.7);
plot(Month_labels,QQG_SM.PercSmelt./100,'Color',[0.6 0.6 0.6],'LineStyle',':','LineWidth',1);
plot(Month_labels,QQG_SM.PercSub./100,'Color',[0.9 0.17 0.71],'LineStyle',':','LineWidth',1);
ba5(1).CData = [0.49 0.73 0.91]; ba5(2).CData = [1 1 1]*0.8; ba5(3).CData = [0.95 0.61 0.73];
title('Quisoquipina Glacier (5180 m)');ylim([0 1]);
yyaxis right; plot(Month_labels,QQG_AM.mean_SND,'Color',[0.05 0.6 0.73],'LineStyle','--','LineWidth',1); set(gca,'YColor',[0.05 0.6 0.73]);ylim([0 0.7]);

Abla_names=["Ice melt","Snow melt","Sublimation","Albedo","Fraction snowmelt","Fraction sublimation","Snow depth"];
f1L=legend(Abla_names,'Location','SouthOutside','Orientation','Horizontal','NumColumns',3);
f1P = [0.5 0.02 0.03 0.03];
set(f1L,'Position', f1P,'Units', 'normalized');

%============ Scatter plots against elevation =====================
%mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ
All_A_fluxes_scatter=[table2array(SG_A_fluxes); table2array(AG_A_fluxes); table2array(CG_A_fluxes); table2array(QQG_A_fluxes); table2array(QIC_A_fluxes)];
All_SWnet_flux_sc = reshape(All_A_fluxes_scatter(:,1),[3,5]); %So here each column is a glacier, and the rows are mean all, dry, wet
All_LWnet_flux_sc = reshape(All_A_fluxes_scatter(:,2),[3,5]); 
All_H_flux_sc = reshape(All_A_fluxes_scatter(:,3),[3,5]); 
All_QE_flux_sc = reshape(All_A_fluxes_scatter(:,4),[3,5]); 
All_Qv_flux_sc = reshape(All_A_fluxes_scatter(:,5),[3,5]); 
All_Qfm_flux_sc = reshape(All_A_fluxes_scatter(:,6),[3,5]); 
All_G_flux_sc = reshape(All_A_fluxes_scatter(:,7),[3,5]); 
All_dQ_flux_sc = reshape(All_A_fluxes_scatter(:,8),[3,5]); 

%Scatter plot for ablation

SG_A_fluxes(1,:)=array2table(SG_A_all(:,[40 41 7 19 21 20 6 35])); %All
SG_A_fluxes(2,:)=SG_A_season(1,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); %Dry season
SG_A_fluxes(3,:)=SG_A_season(2,{'mean_SWnet','mean_LWnet','mean_H','mean_QE','mean_Qv','mean_Qfm','mean_G','mean_dQ'}); %Wet season

SG_A_all_sub=SG_A_all(:,4)+SG_A_all(:,5); %EICE+ESN 
SG_A_all_abl=SG_A_all_sub+SG_A_all(:,10)+SG_A_all(:,26); %EICE+ESN+Imelt+Smelt
AG_A_all_sub=AG_A_all(:,4)+AG_A_all(:,5); %EICE+ESN 
AG_A_all_abl=AG_A_all_sub+AG_A_all(:,10)+AG_A_all(:,26); %EICE+ESN+Imelt+Smelt
CG_A_all_sub=CG_A_all(:,4)+CG_A_all(:,5); %EICE+ESN 
CG_A_all_abl=CG_A_all_sub+CG_A_all(:,10)+CG_A_all(:,25); %EICE+ESN+Imelt+Smelt %Note different column number
QQG_A_all_sub=QQG_A_all(:,4)+QQG_A_all(:,5); %EICE+ESN 
QQG_A_all_abl=QQG_A_all_sub+QQG_A_all(:,10)+QQG_A_all(:,26); %EICE+ESN+Imelt+Smelt
QIC_A_all_sub=QIC_A_all(:,4)+QIC_A_all(:,5); %EICE+ESN 
QIC_A_all_abl=QIC_A_all_sub+QIC_A_all(:,10)+QIC_A_all(:,25); %EICE+ESN+Imelt+Smelt %Note different column number

SG_A_season.mean_Sub=SG_A_season.mean_EICE+SG_A_season.mean_ESN;
SG_A_season.mean_Abl=SG_A_season.mean_Imelt+SG_A_season.mean_Smelt+SG_A_season.mean_EICE+SG_A_season.mean_ESN;
AG_A_season.mean_Sub=AG_A_season.mean_EICE+AG_A_season.mean_ESN;
AG_A_season.mean_Abl=AG_A_season.mean_Imelt+AG_A_season.mean_Smelt+AG_A_season.mean_EICE+AG_A_season.mean_ESN;
CG_A_season.mean_Sub=CG_A_season.mean_EICE+CG_A_season.mean_ESN;
CG_A_season.mean_Abl=CG_A_season.mean_Imelt+CG_A_season.mean_Smelt+CG_A_season.mean_EICE+CG_A_season.mean_ESN;
QQG_A_season.mean_Sub=QQG_A_season.mean_EICE+QQG_A_season.mean_ESN;
QQG_A_season.mean_Abl=QQG_A_season.mean_Imelt+QQG_A_season.mean_Smelt+QQG_A_season.mean_EICE+QQG_A_season.mean_ESN;
QIC_A_season.mean_Sub=QIC_A_season.mean_EICE+QIC_A_season.mean_ESN;
QIC_A_season.mean_Abl=QIC_A_season.mean_Imelt+QIC_A_season.mean_Smelt+QIC_A_season.mean_EICE+QIC_A_season.mean_ESN;

All_ablat_seas(1,1:2)=SG_A_all(:,[10 26]); All_ablat_seas(1,3)=SG_A_all_sub;  All_ablat_seas(1,4)=SG_A_all_abl; %SG
All_ablat_seas(2,1:2)=AG_A_all(:,[10 26]); All_ablat_seas(2,3)=AG_A_all_sub;  All_ablat_seas(2,4)=AG_A_all_abl; %AG
All_ablat_seas(3,1:2)=CG_A_all(:,[10 25]); All_ablat_seas(3,3)=CG_A_all_sub;  All_ablat_seas(3,4)=CG_A_all_abl; %CG %Note different column number
All_ablat_seas(4,1:2)=QQG_A_all(:,[10 26]); All_ablat_seas(4,3)=QQG_A_all_sub;  All_ablat_seas(4,4)=QQG_A_all_abl; %QQG
All_ablat_seas(5,1:2)=QIC_A_all(:,[10 25]); All_ablat_seas(5,3)=QIC_A_all_sub;  All_ablat_seas(5,4)=QIC_A_all_abl; %QIC %Note different column number

Dry_ablat_seas(1,:)=SG_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub','mean_Abl'}); %SG
Dry_ablat_seas(2,:)=AG_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub','mean_Abl'}); %SG
Dry_ablat_seas(3,:)=CG_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub','mean_Abl'}); %CG
Dry_ablat_seas(4,:)=QQG_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub','mean_Abl'}); %QQG
Dry_ablat_seas(5,:)=QIC_A_season(1,{'mean_Imelt','mean_Smelt','mean_Sub','mean_Abl'}); %CG

Wet_ablat_seas(1,:)=SG_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub','mean_Abl'}); %SG
Wet_ablat_seas(2,:)=AG_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub','mean_Abl'}); %SG
Wet_ablat_seas(3,:)=CG_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub','mean_Abl'}); %CG
Wet_ablat_seas(4,:)=QQG_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub','mean_Abl'}); %QQG
Wet_ablat_seas(5,:)=QIC_A_season(2,{'mean_Imelt','mean_Smelt','mean_Sub','mean_Abl'}); %CG


%Figure with %snow melt and ice melt per season
%Seasons
SG_S_season.sum_Sub=SG_S_season.sum_ESN+SG_S_season.sum_EICE; %Sublimation
SG_S_season.sum_Abl=SG_S_season.sum_Tmelt+SG_S_season.sum_Sub; %Ablation
SG_S_season.perc_Smelt=(SG_S_season.sum_Smelt./SG_S_season.sum_Abl).*100; %Perc snow melt
SG_S_season.perc_Imelt=(SG_S_season.sum_Imelt./SG_S_season.sum_Abl).*100; %Perc ice melt
SG_S_season.perc_Sub=(SG_S_season.sum_Sub./SG_S_season.sum_Abl).*100; %Perc sublimation
AG_S_season.sum_Sub=AG_S_season.sum_ESN+AG_S_season.sum_EICE; %Sublimation
AG_S_season.sum_Abl=AG_S_season.sum_Tmelt+AG_S_season.sum_Sub; %Ablation
AG_S_season.perc_Smelt=(AG_S_season.sum_Smelt./AG_S_season.sum_Abl).*100; %Perc snow melt
AG_S_season.perc_Imelt=(AG_S_season.sum_Imelt./AG_S_season.sum_Abl).*100; %Perc ice melt
AG_S_season.perc_Sub=(AG_S_season.sum_Sub./AG_S_season.sum_Abl).*100; %Perc sublimation
CG_S_season.sum_Sub=CG_S_season.sum_ESN+CG_S_season.sum_EICE; %Sublimation
CG_S_season.sum_Abl=CG_S_season.sum_Tmelt+CG_S_season.sum_Sub; %Ablation
CG_S_season.perc_Smelt=(CG_S_season.sum_Smelt./CG_S_season.sum_Abl).*100; %Perc snow melt
CG_S_season.perc_Imelt=(CG_S_season.sum_Imelt./CG_S_season.sum_Abl).*100; %Perc ice melt
CG_S_season.perc_Sub=(CG_S_season.sum_Sub./CG_S_season.sum_Abl).*100; %Perc sublimation
QQG_S_season.sum_Sub=QQG_S_season.sum_ESN+QQG_S_season.sum_EICE; %Sublimation
QQG_S_season.sum_Abl=QQG_S_season.sum_Tmelt+QQG_S_season.sum_Sub; %Ablation
QQG_S_season.perc_Smelt=(QQG_S_season.sum_Smelt./QQG_S_season.sum_Abl).*100; %Perc snow melt
QQG_S_season.perc_Imelt=(QQG_S_season.sum_Imelt./QQG_S_season.sum_Abl).*100; %Perc ice melt
QQG_S_season.perc_Sub=(QQG_S_season.sum_Sub./QQG_S_season.sum_Abl).*100; %Perc sublimation
QIC_S_season.sum_Sub=QIC_S_season.sum_ESN+QIC_S_season.sum_EICE; %Sublimation
QIC_S_season.sum_Abl=QIC_S_season.sum_Tmelt+QIC_S_season.sum_Sub; %Ablation
QIC_S_season.perc_Smelt=(QIC_S_season.sum_Smelt./QIC_S_season.sum_Abl).*100; %Perc snow melt
QIC_S_season.perc_Imelt=(QIC_S_season.sum_Imelt./QIC_S_season.sum_Abl).*100; %Perc ice melt
QIC_S_season.perc_Sub=(QIC_S_season.sum_Sub./QIC_S_season.sum_Abl).*100; %Perc sublimation

Dry_ablat_seas_perc(1,:)=SG_S_season(1,{'perc_Imelt','perc_Smelt','perc_Sub'}); %SG
Dry_ablat_seas_perc(2,:)=AG_S_season(1,{'perc_Imelt','perc_Smelt','perc_Sub'}); %AG
Dry_ablat_seas_perc(3,:)=CG_S_season(1,{'perc_Imelt','perc_Smelt','perc_Sub'}); %CG
Dry_ablat_seas_perc(4,:)=QQG_S_season(1,{'perc_Imelt','perc_Smelt','perc_Sub'}); %QQG
Dry_ablat_seas_perc(5,:)=QIC_S_season(1,{'perc_Imelt','perc_Smelt','perc_Sub'}); %QIC

Wet_ablat_seas_perc(1,:)=SG_S_season(2,{'perc_Imelt','perc_Smelt','perc_Sub'}); %SG
Wet_ablat_seas_perc(2,:)=AG_S_season(2,{'perc_Imelt','perc_Smelt','perc_Sub'}); %AG
Wet_ablat_seas_perc(3,:)=CG_S_season(2,{'perc_Imelt','perc_Smelt','perc_Sub'}); %CG
Wet_ablat_seas_perc(4,:)=QQG_S_season(2,{'perc_Imelt','perc_Smelt','perc_Sub'}); %QQG
Wet_ablat_seas_perc(5,:)=QIC_S_season(2,{'perc_Imelt','perc_Smelt','perc_Sub'}); %QIC

%All
SG_T_labels=SG_T.Properties.VariableNames(2:end);
SG_S_all_T=array2table(SG_S_all,'VariableNames',SG_T_labels);%Creating a table
SG_S_all_T.Sub=SG_S_all_T.ESN+SG_S_all_T.EICE; %Sublimation
SG_S_all_T.Abl=SG_S_all_T.Tmelt+SG_S_all_T.Sub; %Ablation
SG_S_all_T.perc_Smelt=(SG_S_all_T.Smelt./SG_S_all_T.Abl).*100; %Perc snow melt
SG_S_all_T.perc_Imelt=(SG_S_all_T.Imelt./SG_S_all_T.Abl).*100; %Perc ice melt
SG_S_all_T.perc_Sub=(SG_S_all_T.Sub./SG_S_all_T.Abl).*100; %Perc sublimation
AG_T_labels=AG_T.Properties.VariableNames(2:end);
AG_S_all_T=array2table(AG_S_all,'VariableNames',AG_T_labels);%Creating a table
AG_S_all_T.Sub=AG_S_all_T.ESN+AG_S_all_T.EICE; %Sublimation
AG_S_all_T.Abl=AG_S_all_T.Tmelt+AG_S_all_T.Sub; %Ablation
AG_S_all_T.perc_Smelt=(AG_S_all_T.Smelt./AG_S_all_T.Abl).*100; %Perc snow melt
AG_S_all_T.perc_Imelt=(AG_S_all_T.Imelt./AG_S_all_T.Abl).*100; %Perc ice melt
AG_S_all_T.perc_Sub=(AG_S_all_T.Sub./AG_S_all_T.Abl).*100; %Perc sublimation
CG_T_labels=CG_T.Properties.VariableNames(2:end);
CG_S_all_T=array2table(CG_S_all,'VariableNames',CG_T_labels);%Creating a table
CG_S_all_T.Sub=CG_S_all_T.ESN+CG_S_all_T.EICE; %Sublimation
CG_S_all_T.Abl=CG_S_all_T.Tmelt+CG_S_all_T.Sub; %Ablation
CG_S_all_T.perc_Smelt=(CG_S_all_T.Smelt./CG_S_all_T.Abl).*100; %Perc snow melt
CG_S_all_T.perc_Imelt=(CG_S_all_T.Imelt./CG_S_all_T.Abl).*100; %Perc ice melt
CG_S_all_T.perc_Sub=(CG_S_all_T.Sub./CG_S_all_T.Abl).*100; %Perc sublimation
QQG_T_labels=QQG_T.Properties.VariableNames(2:end);
QQG_S_all_T=array2table(QQG_S_all,'VariableNames',QQG_T_labels);%Creating a table
QQG_S_all_T.Sub=QQG_S_all_T.ESN+QQG_S_all_T.EICE; %Sublimation
QQG_S_all_T.Abl=QQG_S_all_T.Tmelt+QQG_S_all_T.Sub; %Ablation
QQG_S_all_T.perc_Smelt=(QQG_S_all_T.Smelt./QQG_S_all_T.Abl).*100; %Perc snow melt
QQG_S_all_T.perc_Imelt=(QQG_S_all_T.Imelt./QQG_S_all_T.Abl).*100; %Perc ice melt
QQG_S_all_T.perc_Sub=(QQG_S_all_T.Sub./QQG_S_all_T.Abl).*100; %Perc sublimation
QIC_T_labels=QIC_T.Properties.VariableNames(2:end);
QIC_S_all_T=array2table(QIC_S_all,'VariableNames',QIC_T_labels);%Creating a table
QIC_S_all_T.Sub=QIC_S_all_T.ESN+QIC_S_all_T.EICE; %Sublimation
QIC_S_all_T.Abl=QIC_S_all_T.Tmelt+QIC_S_all_T.Sub; %Ablation
QIC_S_all_T.perc_Smelt=(QIC_S_all_T.Smelt./QIC_S_all_T.Abl).*100; %Perc snow melt
QIC_S_all_T.perc_Imelt=(QIC_S_all_T.Imelt./QIC_S_all_T.Abl).*100; %Perc ice melt
QIC_S_all_T.perc_Sub=(QIC_S_all_T.Sub./QIC_S_all_T.Abl).*100; %Perc sublimation

All_ablat_seas_perc(1,:)=SG_S_all_T(1,{'perc_Imelt','perc_Smelt','perc_Sub'}); %SG
All_ablat_seas_perc(2,:)=AG_S_all_T(1,{'perc_Imelt','perc_Smelt','perc_Sub'}); %AG
All_ablat_seas_perc(3,:)=CG_S_all_T(1,{'perc_Imelt','perc_Smelt','perc_Sub'}); %CG
All_ablat_seas_perc(4,:)=QQG_S_all_T(1,{'perc_Imelt','perc_Smelt','perc_Sub'}); %QQG
All_ablat_seas_perc(5,:)=QIC_S_all_T(1,{'perc_Imelt','perc_Smelt','perc_Sub'}); %QIC

%Combined figure
Ablat_colours(1,:)=[0.49 0.73 0.91]; %ice 
Ablat_colours(2,:)=[0.6 0.6 0.6]; %snow
Ablat_colours(3,:)=[0.95 0.61 0.73]; %sublimation
Ablat_colours(4,:)=[0 0 0]; %ablation

%Combined figure - tight
figure20=figure();
f=tiledlayout(2,3,'TileSpacing','tight');
ax=nexttile; ax.ColorOrder = Ablat_colours; hold on; 
plot(All_el,All_ablat_seas(:,1),'-d',All_el,All_ablat_seas(:,2),'-d',All_el,All_ablat_seas(:,3),'-d',All_el,All_ablat_seas(:,4),'-d');
ylabel('Mean ablation (mm h^-^1)','FontSize',12); title('All'); ylim([-0.1 1.1]); xlim([4700 5700]);
ax=nexttile;  ax.ColorOrder = Ablat_colours; hold on; 
plot(All_el,Dry_ablat_seas.mean_Imelt,'-d',All_el,Dry_ablat_seas.mean_Smelt,'-d',All_el,Dry_ablat_seas.mean_Sub,'-d',All_el,Dry_ablat_seas.mean_Abl,'-d');
title('Dry season'); ylim([-0.1 1.1]); xlim([4700 5700]);
ax=nexttile; ax.ColorOrder = Ablat_colours; hold on; 
plot(All_el,Wet_ablat_seas.mean_Imelt,'-d',All_el,Wet_ablat_seas.mean_Smelt,'-d',All_el,Wet_ablat_seas.mean_Sub,'-d',All_el,Wet_ablat_seas.mean_Abl,'-d');
title('Wet season'); ylim([-0.1 1.1]); xlim([4700 5700]);
legend('Ice melt','Snow melt','Sublimation','Ablation');
%Percentages
ax=nexttile; ax.ColorOrder = Ablat_colours; hold on;  
plot(All_el,All_ablat_seas_perc.perc_Imelt,'-d',All_el,All_ablat_seas_perc.perc_Smelt,'-d',All_el,All_ablat_seas_perc.perc_Sub,'-d');
ylabel('Ablation (%)','FontSize',12); title('All'); ylim([0 100]);xlim([4700 5700]);
yyaxis right
plot(All_el,All_ablat_seas(:,4),'--d');
ax=nexttile; ax.ColorOrder = Ablat_colours; hold on; 
plot(All_el,Dry_ablat_seas_perc.perc_Imelt,'-d',All_el,Dry_ablat_seas_perc.perc_Smelt,'-d',All_el,Dry_ablat_seas_perc.perc_Sub,'-d');
xlabel('Elevation (m a.s.l.)','FontSize',12); title('Dry season'); ylim([0 100]);xlim([4700 5700]);
yyaxis right
plot(All_el,Dry_ablat_seas.mean_Abl,'--d');
ax=nexttile; ax.ColorOrder = Ablat_colours; hold on; 
plot(All_el,Wet_ablat_seas_perc.perc_Imelt,'-d',All_el,Wet_ablat_seas_perc.perc_Smelt,'-d',All_el,Wet_ablat_seas_perc.perc_Sub,'-d');
title('Wet season'); ylim([0 100]); xlim([4700 5700]);
yyaxis right
plot(All_el,Wet_ablat_seas.mean_Abl,'--d');ylabel('Mean ablation (mm w.e. h^-^1)','FontSize',12);
legend('Ice melt (%)','Snow melt (%)','Sublimation (%)','Mean ablation');
% Create textbox
annotation(figure20,'textbox',[0.07 0.921153846153847 0.0184621941978239 0.0465625456137788],'String',{'a'},'LineStyle','none','FontSize',14,'FitBoxToText','off');
annotation(figure20,'textbox',[0.07 0.465384615384617 0.0184621941978239 0.0465625456137788],'String','b','LineStyle','none','FontSize',14,'FitBoxToText','off');
annotation(figure20,'textbox',[0.577514678326046 0.622217578840234 0.0352273579013117 0.0359937402190923],'String','QIC','FitBoxToText','off','EdgeColor','none');
annotation(figure20,'textbox',[0.500687695190506 0.683163849890755 0.0352273579013117 0.0359937402190923],'String','QQG','FitBoxToText','off','EdgeColor','none');
annotation(figure20,'textbox',[0.436352904434728 0.744110120941276 0.0352273579013118 0.0359937402190923],'String','CG','FitBoxToText','off','EdgeColor','none');
annotation(figure20,'textbox',[0.432605246720799 0.81788718589717 0.0352273579013119 0.0359937402190923],'String','AG','FitBoxToText','off','EdgeColor','none');
annotation(figure20,'textbox',[0.431980637101811 0.858785341470547 0.0352273579013118 0.0359937402190923],'String',{'SG'},'FitBoxToText','off','EdgeColor','none');

%====== El nino/ La nina analysis ======================

%Need to calculate 3 monthly running mean
%Will use monthly means as a starting point

%Load NOAA ONI data
load('NOAA_ONI.mat','NOAA_ONI');
NOAAt1=datetime(1950,01,01,0,0,0);
NOAAt2=datetime(2020,07,01,0,0,0);
NOAA_ONI.months=(NOAAt1:calmonths(1):NOAAt2)'; %Note the month is the middle month of the 3 month series
NOAA_ONI.numonths=datenum(NOAA_ONI.months); %Can use month number to join to data
NOAA_size=size(NOAA_ONI,1);

%===== Shallap Glacier ==============
%Average data
SG_YearMonth=datetime(SG_T.Year,SG_T.Month,1,0,0,0);%This is the year and  month but no day and time (given first day so doesn't give month before)
SG_T.YearMonth=datenum(SG_YearMonth); %As a number
%Now calculate average per month based on each month
SG_ApM=grpstats(SG_T(:,2:end),'YearMonth','mean'); %Calculates mean per month (for each month, not all together)
SG_Amm3M=movmean(SG_ApM{:,:},3,'omitnan'); %3 month moving average %Good that works, but does wierd things to month and year
SG_Amm3M=[SG_ApM.YearMonth, SG_Amm3M]; %Add correct months back on
SG_3Msize=size(SG_Amm3M,1);
SG_NOAA_ONI=zeros(1,SG_3Msize);
for x=1:NOAA_size
    for y=1:SG_3Msize
        if NOAA_ONI.numonths(x)==SG_Amm3M(y,1)
            SG_NOAA_ONI(y)=NOAA_ONI.ANOM(x);
        end
    end
end

SG_NOAA_ONI=SG_NOAA_ONI'; %Nice that works
%Turn into table and add on NOAA_ONI
SG_table_names=SG_ApM.Properties.VariableNames;
SG_Amm3M_T=array2table(SG_Amm3M,'VariableNames',[{'Correct_Year_Month'} SG_table_names]);
SG_Amm3M_T.NOAA_ONI=SG_NOAA_ONI;
SG_Amm3M_T.mean_Month=SG_ApM.mean_Month; %Corrects months for sorting later
SG_Amm3M_T.mean_Year=SG_ApM.mean_Year; %Corrects year for sorting later
%Remove problem rows
SG_Amm3M_T([1,17,18,end],:) = []; %Note this removes the first and last row as incomplete and the 2 months over the split for SG only (Nov 2011 and May 2012) 
%Note the new edge rows may not have a full 3 months, but will have at
%least 2 months of data
%Split by month
for m=1:12
SG_Amm3M_T_split(m).months=SG_Amm3M_T(SG_Amm3M_T.mean_Month == (m),:); %NICE! NOTE that the index goes BEFORE the dot and splits into 12 structures
end

%===== Artesonraju Glacier ==============
%Average data
AG_YearMonth=datetime(AG_T.Year,AG_T.Month,1,0,0,0);%This is the year and  month but no day and time (given first day so doesn't give month before)
AG_T.YearMonth=datenum(AG_YearMonth); %As a number
%Now calculate average per month based on each month
AG_ApM=grpstats(AG_T(:,2:end),'YearMonth','mean'); %Calculates mean per month (for each month, not all together)
AG_Amm3M=movmean(AG_ApM{:,:},3,'omitnan'); %3 month moving average %Good that works, but does wierd things to month and year
AG_Amm3M=[AG_ApM.YearMonth, AG_Amm3M]; %Add correct months back on
AG_3Msize=size(AG_Amm3M,1);
AG_NOAA_ONI=zeros(1,AG_3Msize);
for x=1:NOAA_size
    for y=1:AG_3Msize
        if NOAA_ONI.numonths(x)==AG_Amm3M(y,1)
            AG_NOAA_ONI(y)=NOAA_ONI.ANOM(x);
        end
    end
end

AG_NOAA_ONI=AG_NOAA_ONI'; %Nice that works
%Turn into table and add on NOAA_ONI
AG_table_names=AG_ApM.Properties.VariableNames;
AG_Amm3M_T=array2table(AG_Amm3M,'VariableNames',[{'Correct_Year_Month'} AG_table_names]);
AG_Amm3M_T.NOAA_ONI=AG_NOAA_ONI;AG_Amm3M_T.mean_Year=AG_ApM.mean_Year; %Corrects year for sorting later
AG_Amm3M_T.mean_Month=AG_ApM.mean_Month; %Corrects months for sorting later

%Remove problem rows
AG_Amm3M_T([1,end],:) = []; %Note this removes the first and last row as incomplete 
%Note the new edge rows may not have a full 3 months, but will have at
%least 2 months of data
%Split by month
for m=1:12
AG_Amm3M_T_split(m).months=AG_Amm3M_T(AG_Amm3M_T.mean_Month == (m),:); %NICE! NOTE that the index goes BEFORE the dot and splits into 12 structures
end

%===== Cuchillacocha Glacier ==============
%Average data
CG_YearMonth=datetime(CG_T.Year,CG_T.Month,1,0,0,0);%This is the year and  month but no day and time (given first day so doesn't give month before)
CG_T.YearMonth=datenum(CG_YearMonth); %As a number
%Now calculate average per month based on each month
CG_ApM=grpstats(CG_T(:,2:end),'YearMonth','mean'); %Calculates mean per month (for each month, not all together)
CG_Amm3M=movmean(CG_ApM{:,:},3,'omitnan'); %3 month moving average %Good that works, but does wierd things to month and year
CG_Amm3M=[CG_ApM.YearMonth, CG_Amm3M]; %Add correct months back on
CG_3Msize=size(CG_Amm3M,1);
CG_NOAA_ONI=zeros(1,CG_3Msize);
for x=1:NOAA_size
    for y=1:CG_3Msize
        if NOAA_ONI.numonths(x)==CG_Amm3M(y,1)
            CG_NOAA_ONI(y)=NOAA_ONI.ANOM(x);
        end
    end
end

CG_NOAA_ONI=CG_NOAA_ONI'; %Nice that works
%Turn into table and add on NOAA_ONI
CG_table_names=CG_ApM.Properties.VariableNames;
CG_Amm3M_T=array2table(CG_Amm3M,'VariableNames',[{'Correct_Year_Month'} CG_table_names]);
CG_Amm3M_T.NOAA_ONI=CG_NOAA_ONI;
CG_Amm3M_T.mean_Month=CG_ApM.mean_Month; %Corrects months for sorting later
%Remove problem rows
CG_Amm3M_T([1,end],:) = []; %Note this removes the first and last row as incomplete 
%Note the new edge rows may not have a full 3 months, but will have at
%least 2 months of data
%Split by month
for m=1:12
CG_Amm3M_T_split(m).months=CG_Amm3M_T(CG_Amm3M_T.mean_Month == (m),:); %NICE! NOTE that the index goes BEFORE the dot and splits into 12 structures
end

%===== Quelccaya Ice Cap ==============
%Average data
QIC_YearMonth=datetime(QIC_T.Year,QIC_T.Month,1,0,0,0);%This is the year and  month but no day and time (given first day so doesn't give month before)
QIC_T.YearMonth=datenum(QIC_YearMonth); %As a number
%Now calculate average per month based on each month
QIC_ApM=grpstats(QIC_T(:,2:end),'YearMonth','mean'); %Calculates mean per month (for each month, not all together)
QIC_Amm3M=movmean(QIC_ApM{:,:},3,'omitnan'); %3 month moving average %Good that works, but does wierd things to month and year
QIC_Amm3M=[QIC_ApM.YearMonth, QIC_Amm3M]; %Add correct months back on
QIC_3Msize=size(QIC_Amm3M,1);
QIC_NOAA_ONI=zeros(1,QIC_3Msize);
for x=1:NOAA_size
    for y=1:QIC_3Msize
        if NOAA_ONI.numonths(x)==QIC_Amm3M(y,1)
            QIC_NOAA_ONI(y)=NOAA_ONI.ANOM(x);
        end
    end
end

QIC_NOAA_ONI=QIC_NOAA_ONI'; %Nice that works
%Turn into table and add on NOAA_ONI
QIC_table_names=QIC_ApM.Properties.VariableNames;
QIC_Amm3M_T=array2table(QIC_Amm3M,'VariableNames',[{'Correct_Year_Month'} QIC_table_names]);
QIC_Amm3M_T.NOAA_ONI=QIC_NOAA_ONI;
QIC_Amm3M_T.mean_Month=QIC_ApM.mean_Month; %Corrects months for sorting later
QIC_Amm3M_T.mean_Year=QIC_ApM.mean_Year; %Corrects year for sorting later
%Remove problem rows
QIC_Amm3M_T([1,end],:) = []; %Note this removes the first and last row as incomplete 
%Note the new edge rows may not have a full 3 months, but will have at
%least 2 months of data
%Split by month
for m=1:12
QIC_Amm3M_T_split(m).months=QIC_Amm3M_T(QIC_Amm3M_T.mean_Month == (m),:); %NICE! NOTE that the index goes BEFORE the dot and splits into 12 structures
end

%===== Quisoquipina Glacier ==============
%Average data
QQG_YearMonth=datetime(QQG_T.Year,QQG_T.Month,1,0,0,0);%This is the year and  month but no day and time (given first day so doesn't give month before)
QQG_T.YearMonth=datenum(QQG_YearMonth); %As a number
%Now calculate average per month based on each month
QQG_ApM=grpstats(QQG_T(:,2:end),'YearMonth','mean'); %Calculates mean per month (for each month, not all together)
QQG_Amm3M=movmean(QQG_ApM{:,:},3,'omitnan'); %3 month moving average %Good that works, but does wierd things to month and year
QQG_Amm3M=[QQG_ApM.YearMonth, QQG_Amm3M]; %Add correct months back on
QQG_3Msize=size(QQG_Amm3M,1);
QQG_NOAA_ONI=zeros(1,QQG_3Msize);
for x=1:NOAA_size
    for y=1:QQG_3Msize
        if NOAA_ONI.numonths(x)==QQG_Amm3M(y,1)
            QQG_NOAA_ONI(y)=NOAA_ONI.ANOM(x);
        end
    end
end

QQG_NOAA_ONI=QQG_NOAA_ONI'; %Nice that works
%Turn into table and add on NOAA_ONI
QQG_table_names=QQG_ApM.Properties.VariableNames;
QQG_Amm3M_T=array2table(QQG_Amm3M,'VariableNames',[{'Correct_Year_Month'} QQG_table_names]);
QQG_Amm3M_T.NOAA_ONI=QQG_NOAA_ONI;
QQG_Amm3M_T.mean_Month=QQG_ApM.mean_Month; %Corrects months for sorting later
QQG_Amm3M_T.mean_Year=QQG_ApM.mean_Year; %Corrects year for sorting later
%Remove problem rows
QQG_Amm3M_T([1,end],:) = []; %Note this removes the first and last row as incomplete 
%Note the new edge rows may not have a full 3 months, but will have at
%least 2 months of data
%Split by month
for m=1:12
QQG_Amm3M_T_split(m).months=QQG_Amm3M_T(QQG_Amm3M_T.mean_Month == (m),:); %NICE! NOTE that the index goes BEFORE the dot and splits into 12 structures
end



