% All topography QAR plots

clear

max_dist = 27000;
min_t = 7200; % plots start at 2 hour
max_t = 96 * 3600; % hours * seconds in an hour - gives seconds
peak_index = [180; 300; 180; 300; 180; 300];

% Load cross-section ids
id_floc = '/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/';
load(strcat(id_floc,'xs_xs_ids.mat'));
load(strcat(id_floc,'xs_rcl_ids.mat'));

%% Channel down (FP up)
% Collect baseflow maximum for Channel down (FP up)

base30cms_chdown_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chdown/run10001_30cms/'...
    'abq10a12downf30cmsrst01_0000_20160104_100000_his.nc'];
base150cms_chdown_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chdown/run10002_150cms/'...
    'abq10a12downf150cmsrst01_0000_20160103_050000_his.nc'];

% Extract flow at all cross-section times
base30cms_chdown_all_Q = ncread(...
    base30cms_chdown_file,'cross_section_discharge');
base150cms_chdown_all_Q = ncread(...
    base150cms_chdown_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms_chdown_Q = base30cms_chdown_all_Q(:,end);
base150cms_chdown_Q = base150cms_chdown_all_Q(:,end);

clear base30cms_chdown_all_Q base150cms_chdown_all_Q

%% Collect flow under steady peak flow for Channel down (FP up)

peak180cms_chdown_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chdown/run10003_180cms/'...
    'abq10a12chdown180cmsrst01_0000_20160102_230000_his.nc'];
peak300cms_chdown_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chdown/run10004_300cms/'...
    'abq10a12chdown300cmsrst02_0000_20160103_210000_his.nc'];

% Extract flow at all cross-section times
peak180cms_chdown_all_Q = ncread(...
    peak180cms_chdown_file,'cross_section_discharge');
peak300cms_chdown_all_Q = ncread(...
    peak300cms_chdown_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms_chdown_Q = peak180cms_chdown_all_Q(:,end);
peak300cms_chdown_Q = peak300cms_chdown_all_Q(:,end);

clear peak180cms040_all_Q peak300cms040_all_Q

%% Collect peak flow information from Channel down (FP up) flood runs

fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chdown/'];

% Define folder with *his.nc files
fpath = ['run10005_s180p/';'run10006_s300p/';...
    'run10007_m180p/';'run10008_m300p/';...
    'run10009_l180p/';'run10010_l300p/']; % file path

folder_num = numel(fpath(:,1));
xs_chdown_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_chdown_cell(fi,1) = {fpath(fi,:)};
    upstream_peak = peak_index(fi);

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms_chdown_Q;
        xs_base_Q = base150cms_chdown_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms_chdown_Q;
        xs_base_Q = base30cms_chdown_Q;
    end
    
    fpath_fi = strcat(fpath_root,fpath(fi,:));
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

% Extract cross-sectional information

    % Returns cell array structured as the following
    % xs_cell
    % Column 1 - name of model run (assigned earlier)
    % Column 2 - Model time
    % Column 3 - XS x and y coordinates
    % Column 4 - All XS measured flow at all times
    % Column 5 - XS #, Time of max discharge, Maximum discharge
    % Column 6 - XS distance
    % Column 7 - Fraction of upstream peak discharge
    % Column 8 - Flow above baseflow
    
    [xs_chdown_cell(fi,2:10), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his_chap2(fnum, fnames, fpath_fi, max_dist, min_t,...
        max_t, xs_peak_Q, xs_base_Q, upstream_peak);
    
end

%% Channel up (FP down)
% Collect baseflow maximum for Channel up (FP down)

base30cms_chup_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chup/run11001_30cms/'...
    'abq10a12upf30cmsrst01_0000_20160104_040000_his.nc'];
base150cms_chup_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chup/run11002_150cms/'...
    'abq10a12upf150cmsrst02_0000_20160104_190000_his.nc'];

% Extract flow at all cross-section times
base30cms_chup_all_Q = ncread(...
    base30cms_chup_file,'cross_section_discharge');
base150cms_chup_all_Q = ncread(...
    base150cms_chup_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms_chup_Q = base30cms_chup_all_Q(:,end);
base150cms_chup_Q = base150cms_chup_all_Q(:,end);

clear base30cms_chup_all_Q base150cms_chup_all_Q

%% Collect flow under steady peak flow for Channel up (FP down)

peak180cms_chup_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chup/run11003_180cms/'...
    'abq10a12upf180cmsrst02_0000_20160104_110000_his.nc'];
peak300cms_chup_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chup/run11004_300cms/'...
    'abq10a12upf300cmsrst02_0000_20160103_210000_his.nc'];

% Extract flow at all cross-section times
peak180cms_chup_all_Q = ncread(...
    peak180cms_chup_file,'cross_section_discharge');
peak300cms_chup_all_Q = ncread(...
    peak300cms_chup_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms_chup_Q = peak180cms_chup_all_Q(:,end);
peak300cms_chup_Q = peak300cms_chup_all_Q(:,end);

clear peak180cms_chup_all_Q peak300cms_chup_all_Q

%% Collect peak flow information from Channel up (FP down) flood runs

fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12chup/'];

% Define folder with *his.nc files
fpath = ['run11005_s180p/';'run11006_s300p/';...
    'run11007_m180p/';'run11008_m300p/';...
    'run11009_l180p/';'run11010_l300p/']; % file path

folder_num = numel(fpath(:,1));
xs_chup_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_chup_cell(fi,1) = {fpath(fi,:)};
    upstream_peak = peak_index(fi);

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms_chup_Q;
        xs_base_Q = base150cms_chup_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms_chup_Q;
        xs_base_Q = base30cms_chup_Q;
    end
    
    fpath_fi = strcat(fpath_root,fpath(fi,:));
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

% Extract cross-sectional information
    % Returns cell array structured as the following
    % xs_cell
    % Column 1 - name of model run (assigned earlier)
    % Column 2 - Model time
    % Column 3 - XS x and y coordinates
    % Column 4 - All XS measured flow at all times
    % Column 5 - XS #, Time of max discharge, Maximum discharge
    % Column 6 - XS distance
    % Column 7 - Fraction of upstream peak discharge
    % Column 8 - Flow above baseflow
    
    [xs_chup_cell(fi,2:10), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his_chap2(fnum, fnames, fpath_fi, max_dist, min_t,...
        max_t, xs_peak_Q, xs_base_Q, upstream_peak);
    
end

%% Bank down
% Collect baseflow maximum for Bank down

base30cms_bankdown_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankdown/run12001_30cms/'...
    'abq10a12bankdown30cmsrst01_0000_20160104_080000_his.nc'];
base150cms_bankdown_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankdown/run12002_150cms/'...
    'abq10a12bankdown150cmsrst01_0000_20160103_000000_his.nc'];

% Extract flow at all cross-section times
base30cms_bankdown_all_Q = ncread(...
    base30cms_bankdown_file,'cross_section_discharge');
base150cms_bankdown_all_Q = ncread(...
    base150cms_bankdown_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms_bankdown_Q = base30cms_bankdown_all_Q(:,end);
base150cms_bankdown_Q = base150cms_bankdown_all_Q(:,end);

clear base30cms_bankdown_all_Q base150cms_bankdown_all_Q

%% Collect flow under steady peak flow for Bank down

peak180cms_bankdown_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankdown/run12003_180cms/'...
    'abq10a12bankdownf180rst02_0000_20160104_120000_his.nc'];
peak300cms_bankdown_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankdown/run12004_300cms/'...
    'abq10a12bankdownf300rst02_0000_20160103_200000_his.nc'];

% Extract flow at all cross-section times
peak180cms_bankdown_all_Q = ncread(...
    peak180cms_bankdown_file,'cross_section_discharge');
peak300cms_bankdown_all_Q = ncread(...
    peak300cms_bankdown_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms_bankdown_Q = peak180cms_bankdown_all_Q(:,end);
peak300cms_bankdown_Q = peak300cms_bankdown_all_Q(:,end);

clear peak180cms_bankdown_all_Q peak300cms_bankdown_all_Q

%% Collect peak flow information from Bank down flood runs

fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankdown/'];

% Define folder with *his.nc files
fpath = ['run12005_s180p/';'run12006_s300p/';...
    'run12007_m180p/';'run12008_m300p/';...
    'run12009_l180p/';'run12010_l300p/']; % file path

folder_num = numel(fpath(:,1));
xs_bankdown_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_bankdown_cell(fi,1) = {fpath(fi,:)};
    upstream_peak = peak_index(fi);

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms_bankdown_Q;
        xs_base_Q = base150cms_bankdown_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms_bankdown_Q;
        xs_base_Q = base30cms_bankdown_Q;
    end
    
    fpath_fi = strcat(fpath_root,fpath(fi,:));
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

% Extract cross-sectional information
    % Returns cell array structured as the following
    % xs_cell
    % Column 1 - name of model run (assigned earlier)
    % Column 2 - Model time
    % Column 3 - XS x and y coordinates
    % Column 4 - All XS measured flow at all times
    % Column 5 - XS #, Time of max discharge, Maximum discharge
    % Column 6 - XS distance
    % Column 7 - Fraction of upstream peak discharge
    % Column 8 - Flow above baseflow
    
    [xs_bankdown_cell(fi,2:10), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his_chap2(fnum, fnames, fpath_fi, max_dist, min_t,...
        max_t, xs_peak_Q, xs_base_Q, upstream_peak);
    
end

%% Bank up
% Collect baseflow maximum for Bank up

base30cms_bankup_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankup/run13001_30cms/'...
    'abq10a12bankup30cmsrst01_0000_20160104_030000_his.nc'];
base150cms_bankup_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankup/run13002_150cms/'...
    'abq10a12bankup150cmsrst02_0000_20160104_040000_his.nc'];

% Extract flow at all cross-section times
base30cms_bankup_all_Q = ncread(...
    base30cms_bankup_file,'cross_section_discharge');
base150cms_bankup_all_Q = ncread(...
    base150cms_bankup_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms_bankup_Q = base30cms_bankup_all_Q(:,end);
base150cms_bankup_Q = base150cms_bankup_all_Q(:,end);

clear base30cms_bankup_all_Q base150cms_bankup_all_Q

%% Collect flow under steady peak flow for Bank up

peak180cms_bankup_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankup/run13003_180cms/'...
    'abq10a12bankupf180rst02_0000_20160104_120000_his.nc'];
peak300cms_bankup_file = [...
    '/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankup/run13004_300cms/'...
    'abq10a12bankupf300rst02_0000_20160103_210000_his.nc'];

% Extract flow at all cross-section times
peak180cms_bankup_all_Q = ncread(...
    peak180cms_bankup_file,'cross_section_discharge');
peak300cms_bankup_all_Q = ncread(...
    peak300cms_bankup_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms_bankup_Q = peak180cms_bankup_all_Q(:,end);
peak300cms_bankup_Q = peak300cms_bankup_all_Q(:,end);

clear peak180cms_bankup_all_Q peak300cms_bankup_all_Q

%% Collect peak flow information from Bank up flood runs

fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12bankup/'];

% Define folder with *his.nc files
fpath = ['run13005_s180p/';'run13006_s300p/';...
    'run13007_m180p/';'run13008_m300p/';...
    'run13009_l180p/';'run13010_l300p/']; % file path

folder_num = numel(fpath(:,1));
xs_bankup_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_bankup_cell(fi,1) = {fpath(fi,:)};
    upstream_peak = peak_index(fi);

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms_bankup_Q;
        xs_base_Q = base150cms_bankup_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms_bankup_Q;
        xs_base_Q = base30cms_bankup_Q;
    end
    
    fpath_fi = strcat(fpath_root,fpath(fi,:));
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

% Extract cross-sectional information
    % Returns cell array structured as the following
    % xs_cell
    % Column 1 - name of model run (assigned earlier)
    % Column 2 - Model time
    % Column 3 - XS x and y coordinates
    % Column 4 - All XS measured flow at all times
    % Column 5 - XS #, Time of max discharge, Maximum discharge
    % Column 6 - XS distance
    % Column 7 - Fraction of upstream peak discharge
    % Column 8 - Flow above baseflow
    
    [xs_bankup_cell(fi,2:10), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his_chap2(fnum, fnames, fpath_fi, max_dist, min_t,...
        max_t, xs_peak_Q, xs_base_Q, upstream_peak);
    
end

%% Baseline Conditions
% Collect baseflow maximum for bank/fp n = 0.060

base30cms060_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3047_30cmsbank060/'...
    'abq10a12bf30cmsfp060rst01_0000_20160103_220000_his.nc'];
base150cms060_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3048_150cmsbank060/'...
    'abq10a12bf150cmsfp060rst01_0000_20160103_190000_his.nc'];

% Extract flow at all cross-section times
base30cms060_all_Q = ncread(base30cms060_file,'cross_section_discharge');
base150cms060_all_Q = ncread(base150cms060_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms060_Q = base30cms060_all_Q(1:max_dist/250*3,end);
base150cms060_Q = base150cms060_all_Q(1:max_dist/250*3,end);

clear base30cms060_all_Q base150cms060_all_Q

%% Collect flow under steady peak flow for fp/bank n = 0.060

peak180cms060_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3049_180cmsbank060/'...
    'abq10a12b180cmsrst01_0000_20160102_190000_his.nc'];
peak300cms060_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3050_300cmsbank060/'...
    'abq10a12bf300cmsfp060rst02_0000_20160103_200000_his.nc'];

% Extract flow at all cross-section times
peak180cms060_all_Q = ncread(peak180cms060_file,'cross_section_discharge');
peak300cms060_all_Q = ncread(peak300cms060_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms060_Q = peak180cms060_all_Q(1:max_dist/250*3,end);
peak300cms060_Q = peak300cms060_all_Q(1:max_dist/250*3,end);

clear peak180cms060_all_Q peak300cms060_all_Q

%% Collect peak flow information from fp/bank n = 0.060 flood runs

fpath_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/'];

% Define folder with *his.nc files
fpath = ['run3053_s180fp060/';'run3054_s300fp060/';...
    'run3055_m180fp060/';'run3056_m300fp060/';...
    'run3057_l180fp060/';'run3058_l300fp060/']; % file path

folder_num = numel(fpath(:,1));
xs_fp060_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_fp060_cell(fi,1) = {fpath(fi,:)};
    upstream_peak = peak_index(fi);

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms060_Q;
        xs_base_Q = base150cms060_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms060_Q;
        xs_base_Q = base30cms060_Q;
    end
    
    fpath_fi = strcat(fpath_root,fpath(fi,:));
    
    fnames = struct2cell(dir(strcat(fpath_fi,'*his.nc')))'; 
    [fnum,~] = size(fnames); % number of files

% Extract cross-sectional information
    % Returns cell array structured as the following
    % xs_cell
    % Column 1 - name of model run (assigned earlier)
    % Column 2 - Model time
    % Column 3 - XS x and y coordinates
    % Column 4 - All XS measured flow at all times
    % Column 5 - XS #, Time of max discharge, Maximum discharge
    % Column 6 - XS distance
    % Column 7 - Fraction of upstream peak discharge
    % Column 8 - Flow above baseflow
    
    [xs_fp060_cell(fi,2:10), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his_chap2(fnum, fnames, fpath_fi, max_dist, min_t,...
        max_t, xs_peak_Q, xs_base_Q, upstream_peak);
    
end



%% Plot attenuation ratio
h = figure(1);
h.Units = 'inches';
%h.Position = [1 1 11 15];
h.PaperSize = [8.5 11];
set(gcf,'color','white');
xs_dist = (250:250:max_dist)/1000;
for i = 1:6
        
    fp060_attenuation = xs_fp060_cell{i,9};
    chdown_attenuation = xs_chdown_cell{i,9};
    chup_attenuation = xs_chup_cell{i,9};
    bankdown_attenuation = xs_bankdown_cell{i,9};
    bankup_attenuation = xs_bankup_cell{i,9};
    
    max_dist_id = max_dist/250;
    
    subplot(3,2,i);
    p1=...
        plot(xs_dist(1:max_dist_id),fp060_attenuation(1:max_dist_id,1),'k-');
    hold on
    p2=...
        plot(xs_dist(1:max_dist_id),chdown_attenuation(1:max_dist_id,1),'--');
    p3=plot(xs_dist(1:max_dist_id),chup_attenuation(1:max_dist_id,1),'--');
    ax = gca;
    ax.ColorOrderIndex = 1;
    p4=plot(xs_dist,bankdown_attenuation(1:max_dist/250,1));
    p5=plot(xs_dist,bankup_attenuation(1:max_dist/250,1));
    if i == 1
        ylim([0.7 1]);
    elseif i == 2
        ylim([0.7 1]);
    elseif i == 3
        ylim([0.7 1]);
        ylabel('Discharge Attenuation Ratio');
    elseif i == 4
        ylim([0.7 1]);
    elseif i == 5 
        ylim([0.7 1]);
        xlabel('Downstream distance (km)');
%         legend([p1 p2 p3],'Baseline','CD','CU',...
%             'Location','southwest')
%         legend('boxoff');
    elseif i == 6
        ylim([0.7 1]);
        xlabel('Downstream distance (km)');
%         legend([p4 p5],'BD','BU',...
%             'Location','southwest')
%         legend('boxoff');
    end
    xlim([0 28]);
    set(gca,'Fontsize',14);
    
    % End percentages

    BD_end(i,1) = bankdown_attenuation(max_dist/250,1);
    BU_end(i,1) = bankup_attenuation(max_dist/250,1);
    CD_end(i,1) = chdown_attenuation(max_dist/250,1);
    CU_end(i,1) = chup_attenuation(max_dist/250,1);
    Baseline_end(i,1) = fp060_attenuation(max_dist/250,1);

    percent_BD(i,1) = (BD_end(i,1) - Baseline_end(i,1)) / ...
        (1 - Baseline_end(i,1));
    percent_BU(i,1) = (BU_end(i,1) - Baseline_end(i,1)) / ...
        (1 - Baseline_end(i,1));
    percent_CD(i,1) = (CD_end(i,1) - Baseline_end(i,1)) / ...
        (1 - Baseline_end(i,1));
    percent_CU(i,1) = (CU_end(i,1) - Baseline_end(i,1)) / ...
        (1 - Baseline_end(i,1));

end

% print('-bestfit',['/Users/Colin/Documents/MATLAB/DFlow/'...
%     'influences_chap2/figures/QARtopography8_5x11'],'-dpdf');

max_topo_bank = max([max(percent_BD) max(percent_BU)]);
max_topo_ch = max([max(percent_CD) max(percent_CU)]);
min_topo_bank = min([min(percent_BD) min(percent_BU)]);
min_topo_ch = min([min(percent_CD) min(percent_CU)]);