data = readtable('C:\Users\Kirsty\Documents\project\Thesis_Appendix_Data\Chapter4\Output_MEANDIR\A12b-MEANDIR_results_HistoricSevernCECcorrected.xlsx');
%% Step 1: Preliminary calcs
% a) Convert flow from m3/s to l/s
data.flow_ls = data.flow.*1000; 

% b) Calculate total major river cations (Ca+Mg+Na) in umol/L
data.totalcations_umolL= data.Ca +data.Na + data.Mg;

% c) Convert model output median, 25 and 75percentile fractions of  (X=Ca, Mg, Na) into umol/L from  silicate, carbonate, evaporite and rain 
%(i.e river conc of cation X * model fraction output)
data.Ca_sil_umolL_median = data.Ca.*data.Casilmed;
data.Ca_sil_umolL_25 = data.Ca.*data.Casil25;
data.Ca_sil_umolL_75= data.Ca.*data.Casil75;
data.Na_sil_umolL_median = data.Na.*data.Nasilmed;
data.Na_sil_umolL_25= data.Na.*data.Nasil25;
data.Na_sil_umolL_75= data.Na.*data.Nasil75;
data.Mg_sil_umolL_median = data.Mg.*data.Mgsilmed;
data.Mg_sil_umolL_25= data.Mg.*data.Mgsil25;
data.Mg_sil_umolL_75= data.Mg.*data.Mgsil75;

data.Ca_carb_umolL_median = data.Ca.*data.Cacarbmed;
data.Ca_carb_umolL_25 = data.Ca.*data.Cacarb25;
data.Ca_carb_umolL_75= data.Ca.*data.Cacarb75;
data.Na_carb_umolL_median = data.Na.*data.Nacarbmed;
data.Na_carb_umolL_25= data.Na.*data.Nacarb25;
data.Na_carb_umolL_75= data.Na.*data.Nacarb75;
data.Mg_carb_umolL_median = data.Mg.*data.Mgcarbmed;
data.Mg_carb_umolL_25= data.Mg.*data.Mgcarb25;
data.Mg_carb_umolL_75= data.Mg.*data.Mgcarb75;

data.Ca_evap_umolL_median = data.Ca.*data.Caevapmed;
data.Ca_evap_umolL_25 = data.Ca.*data.Caevap25;
data.Ca_evap_umolL_75= data.Ca.*data.Caevap75;
data.Na_evap_umolL_median = data.Na.*data.Naevapmed;
data.Na_evap_umolL_25= data.Na.*data.Naevap25;
data.Na_evap_umolL_75= data.Na.*data.Naevap75;
data.Mg_evap_umolL_median = data.Mg.*data.Mgevapmed;
data.Mg_evap_umolL_25= data.Mg.*data.Mgevap25;
data.Mg_evap_umolL_75= data.Mg.*data.Mgevap75;

data.Ca_rain_umolL_median = data.Ca.*data.Carainmed;
data.Ca_rain_umolL_25 = data.Ca.*data.Carain25;
data.Ca_rain_umolL_75= data.Ca.*data.Carain75;
data.Na_rain_umolL_median = data.Na.*data.Narainmed;
data.Na_rain_umolL_25= data.Na.*data.Narain25;
data.Na_rain_umolL_75= data.Na.*data.Narain75;
data.Mg_rain_umolL_median = data.Mg.*data.Mgrainmed;
data.Mg_rain_umolL_25= data.Mg.*data.Mgrain25;
data.Mg_rain_umolL_75= data.Mg.*data.Mgrain75;


%% Step 2: Calculate Xsil, Xcarb, Qsil, Qcarb, etc
%a) Calculate the total cations from silicate, carbonate, evaporite and rain 

data.Total_sil_cations_median_umolL = data.Ca_sil_umolL_median + data.Na_sil_umolL_median + data.Mg_sil_umolL_median;
data.Total_sil_cations_25_umolL = data.Ca_sil_umolL_25 + data.Na_sil_umolL_25 + data.Mg_sil_umolL_25;
data.Total_sil_cations_75_umolL = data.Ca_sil_umolL_75 + data.Na_sil_umolL_75 + data.Mg_sil_umolL_75;

data.Total_carb_cations_median_umolL = data.Ca_carb_umolL_median + data.Na_carb_umolL_median + data.Mg_carb_umolL_median;
data.Total_carb_cations_25_umolL = data.Ca_carb_umolL_25 + data.Na_carb_umolL_25 + data.Mg_carb_umolL_25;
data.Total_carb_cations_75_umolL = data.Ca_carb_umolL_75 + data.Na_carb_umolL_75 + data.Mg_carb_umolL_75;

data.Total_evap_cations_median_umolL = data.Ca_evap_umolL_median + data.Na_evap_umolL_median + data.Mg_evap_umolL_median;
data.Total_evap_cations_25_umolL = data.Ca_evap_umolL_25 + data.Na_evap_umolL_25 + data.Mg_evap_umolL_25;
data.Total_evap_cations_75_umolL = data.Ca_evap_umolL_75 + data.Na_evap_umolL_75 + data.Mg_evap_umolL_75;

data.Total_rain_cations_median_umolL = data.Ca_rain_umolL_median + data.Na_rain_umolL_median + data.Mg_rain_umolL_median;
data.Total_rain_cations_25_umolL = data.Ca_rain_umolL_25 + data.Na_rain_umolL_25 + data.Mg_rain_umolL_25;
data.Total_rain_cations_75_umolL = data.Ca_rain_umolL_75 + data.Na_rain_umolL_75 + data.Mg_rain_umolL_75;

%% b) Calculate Xsil, Xcarb, Xevap and Xrain for the median, 25 and 75 percentile ranges
data.Xsilmedian= data.Total_sil_cations_median_umolL./data.totalcations_umolL.*100;
data.Xsil25= data.Total_sil_cations_25_umolL./data.totalcations_umolL.*100;
data.Xsil75= data.Total_sil_cations_75_umolL./data.totalcations_umolL.*100;

data.Xcarbmedian= data.Total_carb_cations_median_umolL./data.totalcations_umolL.*100;
data.Xcarb25= data.Total_carb_cations_25_umolL./data.totalcations_umolL.*100;
data.Xcarb75= data.Total_carb_cations_75_umolL./data.totalcations_umolL.*100;

data.Xevaomedian= data.Total_evap_cations_median_umolL./data.totalcations_umolL.*100;
data.Xevap25= data.Total_evap_cations_25_umolL./data.totalcations_umolL.*100;
data.Xevap75= data.Total_evap_cations_75_umolL./data.totalcations_umolL.*100;

data.Xrainmedian= data.Total_rain_cations_median_umolL./data.totalcations_umolL.*100;
data.Xrain25= data.Total_rain_cations_25_umolL./data.totalcations_umolL.*100;
data.Xrain75= data.Total_rain_cations_75_umolL./data.totalcations_umolL.*100;
%% c) Convert total cations from each endmember from umol/L to mol/S (i.e. Qsil, Qcarb, Qevap, Qrain)
data.Total_sil_cations_median_mols = data.Total_sil_cations_median_umolL.*data.flow_ls./1000000;
data.Total_sil_cations_25_mols = data.Total_sil_cations_25_umolL.*data.flow_ls./1000000;
data.Total_sil_cations_75_mols = data.Total_sil_cations_75_umolL.*data.flow_ls./1000000;

data.Total_carb_cations_median_mols = data.Total_carb_cations_median_umolL.*data.flow_ls./1000000;
data.Total_carb_cations_25_mols = data.Total_carb_cations_25_umolL.*data.flow_ls./1000000;
data.Total_carb_cations_75_mols = data.Total_carb_cations_75_umolL.*data.flow_ls./1000000;

data.Total_evap_cations_median_mols = data.Total_evap_cations_median_umolL.*data.flow_ls./1000000;
data.Total_evap_cations_25_mols = data.Total_evap_cations_25_umolL.*data.flow_ls./1000000;
data.Total_evap_cations_75_mols = data.Total_evap_cations_75_umolL.*data.flow_ls./1000000;

data.Total_rain_cations_median_mols = data.Total_rain_cations_median_umolL.*data.flow_ls./1000000;
data.Total_rain_cations_25_mols = data.Total_rain_cations_25_umolL.*data.flow_ls./1000000;
data.Total_rain_cations_75_mols = data.Total_rain_cations_75_umolL.*data.flow_ls./1000000;

%% d) Calculate the CO2 consumed due to silicate weathering (mol/s) under the assumption that 1.7 moles of CO2 is removed by 1 Ca and Mg and 2 Na

data.CO2consumption_1_7_sil_median = data.flow_ls.*((data.Ca_sil_umolL_median + data.Mg_sil_umolL_median + data.Na_sil_umolL_median./2)*1.7/1000000);
data.CO2consumption_1_7_sil_25 = data.flow_ls.*((data.Ca_sil_umolL_25 + data.Mg_sil_umolL_25 + data.Na_sil_umolL_25./2)*1.7/1000000);
data.CO2consumption_1_7_sil_75 = data.flow_ls.*((data.Ca_sil_umolL_75 + data.Mg_sil_umolL_75 + data.Na_sil_umolL_75./2)*1.7/1000000);

%% e) Calculate the CO2 consumed due to silicate weathering (mol/s) under the assumption that 1.4 moles of CO2 is removed by 1 Ca and Mg and 2 Na

data.CO2consumption_1_4_sil_median = data.flow_ls.*((data.Ca_sil_umolL_median + data.Mg_sil_umolL_median + data.Na_sil_umolL_median./2)*1.4/1000000);
data.CO2consumption_1_4_sil_25 = data.flow_ls.*((data.Ca_sil_umolL_25 + data.Mg_sil_umolL_25 + data.Na_sil_umolL_25./2)*1.4/1000000);
data.CO2consumption_1_4_sil_75 = data.flow_ls.*((data.Ca_sil_umolL_75 + data.Mg_sil_umolL_75 + data.Na_sil_umolL_75./2)*1.4/1000000);


%% f) Calculate the CO2 due to carbonate weathering (mol/s) 
data.CO2consumption_carb_median = data.flow_ls.*((data.Ca_carb_umolL_median + data.Mg_carb_umolL_median + data.Na_carb_umolL_median./2)/1000000);
data.CO2consumption_carb_25 = data.flow_ls.*((data.Ca_carb_umolL_25 + data.Mg_carb_umolL_25 + data.Na_carb_umolL_25./2)/1000000);
data.CO2consumption_carb_75 = data.flow_ls.*((data.Ca_carb_umolL_75 + data.Mg_carb_umolL_75 + data.Na_carb_umolL_75./2)/1000000);

%% g) Find the monthly averages of Xsil, Xcarb, Qsil, Qcarb, CO2sil, CO2carb
clear av
N = data(:,98:end); % Create table for the data we want an average of
N.date = data.Date; % Add date coloumn to new table
for i = 1:12 % for loop for every month
   x = N(month(N.date) == i,1:end); % Select data for specific month
   if height(x) > 1
   vars = x.Properties.VariableNames; %Get variable names
   vars = vars(1:33);
   av(i,:) = mean(table2array(x(:,1:33)), 'omitnan'); %Convert table to an array and find mean. Add row to new variable
   elseif height(x) == 1
       av(i,:) = table2array(x(:,1:33));
   else 
       c = 0*(ones*1:33);
       av(i,:) = c;
       
   end
end
clear x T i
av = array2table(av); % convert back to table
av.Properties.VariableNames = vars; % add in variable names
clear vars
%% 
av.monthvalue= [31; 28; 31; 30; 31; 30; 31; 31; 30; 31; 30; 31];

%% h) calculating the monthly flux of cations from sils and carb, and CO2 consumed 

av.monthly_Qsil_median = av.Total_sil_cations_median_mols.*av.monthvalue.*60*60*24;
av.monthly_Qsil_25 = av.Total_sil_cations_25_mols.*av.monthvalue.*60*60*24;
av.monthly_Qsil_75 = av.Total_sil_cations_75_mols.*av.monthvalue.*60*60*24;

av.monthly_Qcarb_median = av.Total_carb_cations_median_mols.*av.monthvalue.*60*60*24;
av.monthly_Qcarb_25 = av.Total_carb_cations_25_mols.*av.monthvalue.*60*60*24;
av.monthly_Qcarb_75 = av.Total_carb_cations_75_mols.*av.monthvalue.*60*60*24;

av.monthly_CO2_sil_1_7_median= av.CO2consumption_1_7_sil_median.*av.monthvalue.*60*60*24;
av.monthly_CO2_sil_1_7_25= av.CO2consumption_1_7_sil_25.*av.monthvalue.*60*60*24;
av.monthly_CO2_sil_1_7_75= av.CO2consumption_1_7_sil_75.*av.monthvalue.*60*60*24;

av.monthly_CO2_sil_1_4_median= av.CO2consumption_1_4_sil_median.*av.monthvalue.*60*60*24;
av.monthly_CO2_sil_1_4_25= av.CO2consumption_1_4_sil_25.*av.monthvalue.*60*60*24;
av.monthly_CO2_sil_1_4_75= av.CO2consumption_1_4_sil_75.*av.monthvalue.*60*60*24;

av.monthly_CO2_carb_median= av.CO2consumption_carb_median.*av.monthvalue.*60*60*24;
av.monthly_CO2_carb_25= av.CO2consumption_carb_25.*av.monthvalue.*60*60*24;
av.monthly_CO2_carb_75= av.CO2consumption_carb_75.*av.monthvalue.*60*60*24;

%% Calculating annual fluxes (mol/s) and converting Co2 consumption from mol/s into MtCO2yr-1

annual_Qsil_median_flux_mols= sum(av.monthly_Qsil_median);
annual_Qsil_25_flux_mols= sum(av.monthly_Qsil_25);
annual_Qsil_75_flux_mols= sum(av.monthly_Qsil_75);

annual_Qcarb_median_flux_mols= sum(av.monthly_Qcarb_median);
annual_Qcarb_25_flux_mols= sum(av.monthly_Qcarb_25);
annual_Qcarb_75_flux_mols= sum(av.monthly_Qcarb_75);

annual_CO2consumption_sil_median_1_7_MtCO2yr = sum(av.monthly_CO2_sil_1_7_median)*44.01/1000000/1000000;
annual_CO2consumption_sil_25_1_7_MtCO2yr = sum(av.monthly_CO2_sil_1_7_25)*44.01/1000000/1000000;
annual_CO2consumption_sil_75_1_7_MtCO2yr = sum(av.monthly_CO2_sil_1_7_75)*44.01/1000000/1000000;

annual_CO2consumption_sil_median_1_4_MtCO2yr = sum(av.monthly_CO2_sil_1_4_median)*44.01/1000000/1000000;
annual_CO2consumption_sil_25_1_4_MtCO2yr = sum(av.monthly_CO2_sil_1_4_25)*44.01/1000000/1000000;
annual_CO2consumption_sil_75_1_4_MtCO2yr = sum(av.monthly_CO2_sil_1_4_75)*44.01/1000000/1000000;

annual_CO2consumption_carb_median_MtCO2yr = sum(av.monthly_CO2_carb_median)*44.01/1000000/1000000;
annual_CO2consumption_carb_25_MtCO2yr = sum(av.monthly_CO2_carb_25)*44.01/1000000/1000000;
annual_CO2consumption_carb_75_MtCO2yr = sum(av.monthly_CO2_carb_75)*44.01/1000000/1000000;

Results= table(annual_Qsil_median_flux_mols,annual_Qsil_25_flux_mols,annual_Qsil_75_flux_mols, annual_Qcarb_median_flux_mols,annual_Qcarb_25_flux_mols, annual_Qcarb_75_flux_mols, annual_CO2consumption_sil_median_1_7_MtCO2yr, annual_CO2consumption_sil_25_1_7_MtCO2yr, annual_CO2consumption_sil_75_1_7_MtCO2yr, annual_CO2consumption_sil_median_1_4_MtCO2yr , annual_CO2consumption_sil_25_1_4_MtCO2yr, annual_CO2consumption_sil_75_1_4_MtCO2yr, annual_CO2consumption_carb_median_MtCO2yr, annual_CO2consumption_carb_25_MtCO2yr , annual_CO2consumption_carb_75_MtCO2yr);

%clear annual_CO2consumption_carb_75_MtCO2yr annual_CO2consumption_carb_25_MtCO2yr annual_CO2consumption_carb_median_MtCO2yr
%clear annual_CO2consumption_sil_75_1_4_MtCO2yr annual_CO2consumption_sil_25_1_4_MtCO2yr annual_CO2consumption_sil_median_1_4_MtCO2yr
%clear annual_CO2consumption_sil_75_1_7_MtCO2yr annual_CO2consumption_sil_25_1_7_MtCO2yr annual_CO2consumption_sil_median_1_7_MtCO2yr
%clear annual_Qsil_median_flux_mols annual_Qsil_25_flux_mols annual_Qsil_75_flux_mols
%clear annual_Qcarb_median_flux_mols annual_Qcarb_25_flux_mols annual_Qcarb_75_flux_mols

%% Step 2: calculate how maximums are reduced by Nitric acid 
% a) calculations initial set up
nitratereduction = table();                         %%creates a table to store all info
nitratereduction.NO3agri= data.NO3.*0.7;            %%Using DEFRA's estimate that ~70% of river nitrate is agricultural 
nitratereduction.CWRvsSWR= data.Total_carb_cations_median_umolL./data.Total_sil_cations_median_umolL; %%the carbonate to silicate weathering ratio (HNO3 will be split in the ratio that carbs to sils weather)
nitratereduction.Total = nitratereduction.CWRvsSWR + 1; 
nitratereduction.HNO3tosil= (nitratereduction.NO3agri.*(1./nitratereduction.Total))./2; %%The amount of HNO3 to silicate (dividing by 2 because assuming a 1:2 weathering ratio of Ca/Mg : HNO3)
nitratereduction.HNO3tocarb= (nitratereduction.NO3agri.*(nitratereduction.CWRvsSWR./nitratereduction.Total))./2; %%Amount of HNO3 to carbonate

%% b) finding the total acid needed to weather the cations in solution
nitratereduction.Total_acid_sil_median = data.Ca_sil_umolL_median + data.Mg_sil_umolL_median + data.Na_sil_umolL_median./2;
nitratereduction.Total_acid_carb_median = data.Ca_carb_umolL_median + data.Mg_carb_umolL_median + data.Na_carb_umolL_median./2;

% c) calculating the percentage of the total which is HNO3
nitratereduction.PercentageNitrate_tocarbonate= (nitratereduction.HNO3tocarb./nitratereduction.Total_acid_carb_median)*100;
nitratereduction.PercentageNitrate_tosilicate= (nitratereduction.HNO3tosil./nitratereduction.Total_acid_sil_median)*100;

% d) what is the reduced amount of bicarbonate weathering due to the presence
%of nitric acid
nitratereduction.reducedbicarbonate_sil= nitratereduction.Total_acid_sil_median - nitratereduction.HNO3tosil;
nitratereduction.reducedbicarbonate_carb= nitratereduction.Total_acid_carb_median - nitratereduction.HNO3tocarb;

% e) what is the amount of CO2 consumed by bicarbonate
nitratereduction.CO2_consumed_bicarb_sil= (nitratereduction.reducedbicarbonate_sil.* data.flow_ls)./1000000;
nitratereduction.CO2_consumed_bicarb_carb= (nitratereduction.reducedbicarbonate_carb.* data.flow_ls)./1000000;

% f) What is the amount of CO2 released by HNO3 (for carbonate weathering)
nitratereduction.CO2_released_HNO3_carb= (nitratereduction.HNO3tocarb.* data.flow_ls)./1000000;

% g) What is the overall reduced CO2 consumption by weathering due to HNO3
nitratereduction.Overall_CO2_reduced_nitrate_sil= nitratereduction.CO2_consumed_bicarb_sil.*1.7;
nitratereduction.Overall_CO2_reduced_nitrate_carb=nitratereduction.CO2_consumed_bicarb_carb-  nitratereduction.CO2_released_HNO3_carb;

%% Monthly averages nitrate

N= nitratereduction(:,1:end); % Create table for the data we want an average of

N.date = data.Date; % Add date coloumn to new table
for i = 1:12 % for loop for every month
   A = N(month(N.date) == i,1:end); % Select data for specific month
   if height(A) > 1
   vars = A.Properties.VariableNames; %Get variable names
   vars = vars(1:16);
   averagesNitrate(i,:) = (sum(table2array(A(:,1:16)), 'omitnan'))/length(A.Total_acid_sil_median);
   %mean(table2array(A(:,1:16)), 'omitnan');%Convert table to an array and find mean. Add row to new variable
    elseif height(A) == 1
       averagesNitrate(i,:) = table2array(A(:,1:16));
   else 
       c = 0*(ones*1:16);
       averagesNitrate(i,:) = c;
   end
end
averagesNitrate = array2table(averagesNitrate); % convert back to table
averagesNitrate.Properties.VariableNames = vars; % add in variable names
%clear vars
%clear A N

%% Calculating the changes to the annual CO2 flux due to HNO3

averagesNitrate.averagemonthly_reducedCO2_nitrate_sil=(averagesNitrate.Overall_CO2_reduced_nitrate_sil.*av.monthvalue)*60*60*24;
averagesNitrate.averagemonthly_reducedCO2_nitrate_carb=(averagesNitrate.Overall_CO2_reduced_nitrate_carb.*av.monthvalue)*60*60*24;
%%
annual_reduced_CO2sil_nitrate= sum(averagesNitrate.averagemonthly_reducedCO2_nitrate_sil).*44.01/1000000/1000000;
annual_reduced_CO2carb_nitrate=sum(averagesNitrate.averagemonthly_reducedCO2_nitrate_carb).*44.01/1000000/1000000;
Results2= addvars(Results, annual_reduced_CO2sil_nitrate, annual_reduced_CO2carb_nitrate);

%% Step 3: Calculating the reduction in CO2 consumption due to sulphide weathering

sulphidereduction= table();
sulphidereduction.XexcessSO4median = data.SO4absolute_median- (data.SO4absolute_median.*0.1); %%reducing the excess sulphides for biogenic and volcanic emissions (0.07 and 0.03 respectively)
sulphidereduction.XexcessSO425 = data.SO4absolute_25- (data.SO4absolute_25.*0.1);
sulphidereduction.XexcessSO475 = data.SO4absolute_75- (data.SO4absolute_75.*0.1);

sulphidereduction.XexcessSO4median_tosil= sulphidereduction.XexcessSO4median.*(1./nitratereduction.Total);
sulphidereduction.XexcessSO425_tosil= sulphidereduction.XexcessSO425.*(1./nitratereduction.Total);
sulphidereduction.XexcessSO475_tosil= sulphidereduction.XexcessSO475.*(1./nitratereduction.Total);

sulphidereduction.XexcessSO4median_tocarb= sulphidereduction.XexcessSO4median.*(nitratereduction.CWRvsSWR./nitratereduction.Total);
sulphidereduction.XexcessSO425_tocarb= sulphidereduction.XexcessSO425.*(nitratereduction.CWRvsSWR./nitratereduction.Total);
sulphidereduction.XexcessSO475_tocarb= sulphidereduction.XexcessSO475.*(nitratereduction.CWRvsSWR./nitratereduction.Total);

%% 
%percentage acid that is sulphuric
sulphidereduction.Percent_acid_silicate_sulphuric_median= (sulphidereduction.XexcessSO4median_tosil./nitratereduction.Total_acid_sil_median)*100;
sulphidereduction.Percent_acid_silicate_sulphuric_25= (sulphidereduction.XexcessSO425_tosil./nitratereduction.Total_acid_sil_median)*100;
sulphidereduction.Percent_acid_silicate_sulphuric_75= (sulphidereduction.XexcessSO475_tosil./nitratereduction.Total_acid_sil_median)*100;

sulphidereduction.Percent_acid_carbonate_sulphuric_median= (sulphidereduction.XexcessSO4median_tocarb./nitratereduction.Total_acid_carb_median)*100;
sulphidereduction.Percent_acid_carbonate_sulphuric_25= (sulphidereduction.XexcessSO425_tocarb./nitratereduction.Total_acid_carb_median)*100;
sulphidereduction.Percent_acid_carbonate_sulphuric_75= (sulphidereduction.XexcessSO475_tocarb./nitratereduction.Total_acid_carb_median)*100;

%%
%reduced bicarb due to sulphuric
sulphidereduction.reducedbicarb_sil_median = nitratereduction.Total_acid_sil_median - sulphidereduction.XexcessSO4median_tosil;
sulphidereduction.reducedbicarb_sil_25 = nitratereduction.Total_acid_sil_median - sulphidereduction.XexcessSO425_tosil;
sulphidereduction.reducedbicarb_sil_75 = nitratereduction.Total_acid_sil_median - sulphidereduction.XexcessSO475_tosil;

sulphidereduction.reducedbicarb_carb_median = nitratereduction.Total_acid_carb_median - sulphidereduction.XexcessSO4median_tocarb;
sulphidereduction.reducedbicarb_carb_25 = nitratereduction.Total_acid_carb_median - sulphidereduction.XexcessSO425_tocarb;
sulphidereduction.reducedbicarb_carb_75 = nitratereduction.Total_acid_carb_median - sulphidereduction.XexcessSO475_tocarb;

%% Amount of CO2 consumed by the reduced bicarb
sulphidereduction.CO2consumed_reducedbicarb_sil_median= (sulphidereduction.reducedbicarb_sil_median.*data.flow_ls).*1.7/1000000;
sulphidereduction.CO2consumed_reducedbicarb_sil_25= (sulphidereduction.reducedbicarb_sil_25.*data.flow_ls).*1.7/1000000;
sulphidereduction.CO2consumed_reducedbicarb_sil_75= (sulphidereduction.reducedbicarb_sil_75.*data.flow_ls).*1.7/1000000;

sulphidereduction.CO2consumed_reducedbicarb_carb_median= (sulphidereduction.reducedbicarb_carb_median.*data.flow_ls)./1000000;
sulphidereduction.CO2consumed_reducedbicarb_carb_25= (sulphidereduction.reducedbicarb_carb_25.*data.flow_ls)./1000000;
sulphidereduction.CO2consumed_reducedbicarb_carb_75= (sulphidereduction.reducedbicarb_carb_75.*data.flow_ls)./1000000;

%Released amount of CO2 due to sulphuric (for carbonate)
sulphidereduction.ReleasedCO2_median_carb = (sulphidereduction.XexcessSO4median_tocarb.*data.flow_ls)./1000000;
sulphidereduction.ReleasedCO2_25_carb = (sulphidereduction.XexcessSO425_tocarb.*data.flow_ls)./1000000;
sulphidereduction.ReleasedCO2_75_carb = (sulphidereduction.XexcessSO475_tocarb.*data.flow_ls)./1000000;

%Overall amount of reduction in CO2 due to sulphuric acid
%carbonates
sulphidereduction.OverallCO2reduced_carb_median= sulphidereduction.CO2consumed_reducedbicarb_carb_median -sulphidereduction.ReleasedCO2_median_carb; 
sulphidereduction.OverallCO2reduced_carb_25= sulphidereduction.CO2consumed_reducedbicarb_carb_25 -sulphidereduction.ReleasedCO2_25_carb; 
sulphidereduction.OverallCO2reduced_carb_75= sulphidereduction.CO2consumed_reducedbicarb_carb_75 -sulphidereduction.ReleasedCO2_75_carb; 
%Silicates
sulphidereduction.OverallCO2reduced_sil_median= sulphidereduction.CO2consumed_reducedbicarb_sil_median;
sulphidereduction.OverallCO2reduced_sil_25= sulphidereduction.CO2consumed_reducedbicarb_sil_25;
sulphidereduction.OverallCO2reduced_sil_75= sulphidereduction.CO2consumed_reducedbicarb_sil_75;

%% Monthly averages sulphides

S= sulphidereduction(:,1:end); % Create table for the data we want an average of

S.date = data.Date; % Add date coloumn to new table
for i = 1:12 % for loop for every month
   Q = S(month(S.date) == i,1:end); % Select data for specific month
   if height(Q) > 1
   vars = Q.Properties.VariableNames; %Get variable names
   vars = vars(1:36);
   averagessulphides(i,:) = mean(table2array(Q(:,1:36)), 'omitnan');%Convert table to an array and find mean. Add row to new variable
elseif height(Q) == 1
       averagessulphides(i,:) = table2array(Q(:,1:36));
   else 
       c = 0*(ones*1:36);
       averagessulphides(i,:) = c;
   end
end
averagessulphides = array2table(averagessulphides); % convert back to table
averagessulphides.Properties.VariableNames = vars; % add in variable names
clear vars
clear Q S

%%
XSO4absolute = data.SO4absolute_median;
averageXSO4absolutemedian= mean(XSO4absolute);
averageSO4= mean(data.SO4);
averageXSO4absolute25= mean(data.SO4absolute_25);
averageXSO4absolute75= mean(data.SO4absolute_75);

%% Annual flux reduction calculations for sulphides
averagessulphides.averagemonthly_reducedCO2_sulphide_sil_median= (averagessulphides.OverallCO2reduced_sil_median.*av.monthvalue).*60*60*24;
averagessulphides.averagemonthly_reducedCO2_sulphide_sil_25= (averagessulphides.OverallCO2reduced_sil_25.*av.monthvalue).*60*60*24;
averagessulphides.averagemonthly_reducedCO2_sulphide_sil_75= (averagessulphides.OverallCO2reduced_sil_75.*av.monthvalue).*60*60*24;

averagessulphides.averagemonthly_reducedCO2_sulphide_carb_median= (averagessulphides.OverallCO2reduced_carb_median.*av.monthvalue).*60*60*24;
averagessulphides.averagemonthly_reducedCO2_sulphide_carb_25= (averagessulphides.OverallCO2reduced_carb_25.*av.monthvalue).*60*60*24;
averagessulphides.averagemonthly_reducedCO2_sulphide_carb_75= (averagessulphides.OverallCO2reduced_carb_75.*av.monthvalue).*60*60*24;

%% Calculating the annual CO2 reduced due to sulphuric acid

annual_reduced_CO2sil_sulphide_median= sum(averagessulphides.averagemonthly_reducedCO2_sulphide_sil_median).*44.01/1000000/1000000;
annual_reduced_CO2sil_sulphide_25= sum(averagessulphides.averagemonthly_reducedCO2_sulphide_sil_25).*44.01/1000000/1000000;
annual_reduced_CO2sil_sulphide_75= sum(averagessulphides.averagemonthly_reducedCO2_sulphide_sil_75).*44.01/1000000/1000000;

annual_reduced_CO2carb_sulphide_median=sum(averagessulphides.averagemonthly_reducedCO2_sulphide_carb_median).*44.01/1000000/1000000;
annual_reduced_CO2carb_sulphide_25=sum(averagessulphides.averagemonthly_reducedCO2_sulphide_carb_25).*44.01/1000000/1000000;
annual_reduced_CO2carb_sulphide_75=sum(averagessulphides.averagemonthly_reducedCO2_sulphide_carb_75).*44.01/1000000/1000000;
name= "ElyOuse";
Results3 = addvars(Results2, annual_reduced_CO2sil_sulphide_median, annual_reduced_CO2sil_sulphide_25, annual_reduced_CO2sil_sulphide_75, annual_reduced_CO2carb_sulphide_median,annual_reduced_CO2carb_sulphide_25,annual_reduced_CO2carb_sulphide_75, name, averageXSO4absolutemedian, averageSO4, averageXSO4absolute25, averageXSO4absolute75);
clear annual_reduced_CO2sil_sulphide_median annual_reduced_CO2sil_sulphide_25 annual_reduced_CO2sil_sulphide_75 annual_reduced_CO2carb_sulphide_median annual_reduced_CO2carb_sulphide_25 annual_reduced_CO2carb_sulphide_75


%%
writetable(Results3,'C:\Users\Kirsty\Documents\project\Thesis_Appendix_Data\Chapter4\Output_MEANDIR\Results_Matlab_calcs.xlsx');