% All roughness 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'));

%% FP n = 0.040
% Collect baseflow maximum for FP n = 0.040

base30cms040_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3095_30cmsfp040/'...
    'abq10a12bf30cmsfp040rst01_0000_20160104_040000_his.nc'];
base150cms040_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3096_150cmsfp040/'...
    'abq10a12bf150cmsfp040rst01_0000_20160103_000000_his.nc'];

% Extract flow at all cross-section times
base30cms040_all_Q = ncread(base30cms040_file,'cross_section_discharge');
base150cms040_all_Q = ncread(base150cms040_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms040_Q = base30cms040_all_Q(:,end);
base150cms040_Q = base150cms040_all_Q(:,end);

clear base30cms040_all_Q base150cms040_all_Q

%% Collect flow under steady peak flow for FP n = 0.040

peak180cms040_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3107_f180cmsfp040/'...
    'abq10a12bf180cmsfp040rst02_0000_20160104_130000_his.nc'];
peak300cms040_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3108_f300cmsfp040/'...
    'abq10a12bf300cmsfp040rst02_0000_20160103_230000_his.nc'];

% Extract flow at all cross-section times
peak180cms040_all_Q = ncread(peak180cms040_file,'cross_section_discharge');
peak300cms040_all_Q = ncread(peak300cms040_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms040_Q = peak180cms040_all_Q(:,end);
peak300cms040_Q = peak300cms040_all_Q(:,end);

clear peak180cms040_all_Q peak300cms040_all_Q

%% Collect peak flow information from FP n = 0.040 flood runs

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

% Define folder with *his.nc files
fpath = ['run3101_s180fp040/';'run3102_s300fp040/';...
    'run3103_m180fp040/';'run3104_m300fp040/';...
    'run3105_l180fp040/';'run3106_l300fp040/']; % file path

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

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

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms040_Q;
        xs_base_Q = base150cms040_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms040_Q;
        xs_base_Q = base30cms040_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_fp040_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

%% FP n = 0.060
% Collect baseflow maximum for 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/run3125_150cmsfp060/'...
    'abq10a12bf150fp060rst02_0000_20160104_220000_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(:,end);
base150cms060_Q = base150cms060_all_Q(:,end);

clear base30cms060_all_Q base150cms060_all_Q

%% Collect flow under steady peak flow for FP 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(:,end);
peak300cms060_Q = peak300cms060_all_Q(:,end);

clear peak180cms060_all_Q peak300cms060_all_Q

%% Collect peak flow information from FP n = 0.060 flood runs

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

% Define folder with *his.nc files
fpath = ['run3053_s180fp060/';'run3126_s300fp060/';...
    'run3055_m180fp060/';'run3127_m300fp060/';...
    'run3057_l180fp060/';'run3128_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

%% FP n = 0.080
% Collect baseflow maximum for FP n = 0.080

base30cms080_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3097_30cmsfp080/'...
    'abq10a12bf30cmsfp080rst01_0000_20160104_070000_his.nc'];
base150cms080_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3098_150cmsfp080/'...
    'abq10a12bf150cmsfp080rst01_0000_20160103_000000_his.nc'];

% Extract flow at all cross-section times
base30cms080_all_Q = ncread(base30cms080_file,'cross_section_discharge');
base150cms080_all_Q = ncread(base150cms080_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms080_Q = base30cms080_all_Q(:,end);
base150cms080_Q = base150cms080_all_Q(:,end);

clear base30cms080_all_Q base150cms080_all_Q

%% Collect flow under steady peak flow for FP n = 0.080

peak180cms080_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3109_f180cmsfp080/'...
    'abq10a12bf180cmsfp080rst02_0000_20160104_120000_his.nc'];
peak300cms080_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3110_f300cmsfp080/'...
    'abq10a12bf300cmsfp080rst02_0000_20160103_200000_his.nc'];

% Extract flow at all cross-section times
peak180cms080_all_Q = ncread(peak180cms080_file,'cross_section_discharge');
peak300cms080_all_Q = ncread(peak300cms080_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms080_Q = peak180cms080_all_Q(:,end);
peak300cms080_Q = peak300cms080_all_Q(:,end);

clear peak180cms080_all_Q peak300cms080_all_Q

%% Collect peak flow information from FP n = 0.080 flood runs

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

% Define folder with *his.nc files
fpath = ['run3113_s180fp080/';'run3114_s300fp080/';...
    'run3115_m180fp080/';'run3116_m300fp080/';...
    'run3117_l180fp080/';'run3118_l300fp080/']; % file path

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

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

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms080_Q;
        xs_base_Q = base150cms080_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms080_Q;
        xs_base_Q = base30cms080_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_fp080_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

%% FP n = 0.100
% Collect baseflow maximum for FP n = 0.100

base30cms100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3099_30cmsfp100/'...
    'abq10a12bf30cmsfp100rst01_0000_20160104_100000_his.nc'];
base150cms100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3100_150cmsfp100/'...
    'abq10a12bf150cmsfp100rst02_0000_20160104_210000_his.nc'];

% Extract flow at all cross-section times
base30cms100_all_Q = ncread(base30cms100_file,'cross_section_discharge');
base150cms100_all_Q = ncread(base150cms100_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms100_Q = base30cms100_all_Q(:,end);
base150cms100_Q = base150cms100_all_Q(:,end);

clear base30cms100_all_Q base150cms100_all_Q

%% Collect flow under steady peak flow for FP n = 0.100

peak180cms100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3111_f180cmsfp100/'...
    'abq10a12bf180cmsfp100rst02_0000_20160104_110000_his.nc'];
peak300cms100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3112_f300cmsfp100/'...
    'abq10a12bf300cmsfp100rst03_0000_20160104_150000_his.nc'];

% Extract flow at all cross-section times
peak180cms100_all_Q = ncread(peak180cms100_file,'cross_section_discharge');
peak300cms100_all_Q = ncread(peak300cms100_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms100_Q = peak180cms100_all_Q(:,end);
peak300cms100_Q = peak300cms100_all_Q(:,end);

clear peak180cms080_all_Q peak300cms080_all_Q

%% Collect peak flow information from FP n = 0.100 flood runs

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

% Define folder with *his.nc files
fpath = ['run3119_s180fp100/';'run3120_s300fp100/';...
    'run3121_m180fp100/';'run3122_m300fp100/';...
    'run3123_l180fp100/';'run3124_l300fp100/']; % file path

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

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

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms100_Q;
        xs_base_Q = base150cms100_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms100_Q;
        xs_base_Q = base30cms100_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_fp100_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 n = 0.040
% Collect baseflow maximum for bank n = 0.040

base30cms040_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3065_30cmsbank040/'...
    'abq10a12bf30cmsbank040rst01_0000_20160104_040000_his.nc'];
base150cms040_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3066_150cmsbank040/'...
    'abq10a12bf150cmsbank040rst01_0000_20160103_100000_his.nc'];

% Extract flow at all cross-section times
base30cms040_all_Q = ncread(base30cms040_file,'cross_section_discharge');
base150cms040_all_Q = ncread(base150cms040_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms040_Q = base30cms040_all_Q(1:max_dist/250*3,end);
base150cms040_Q = base150cms040_all_Q(1:max_dist/250*3,end);

clear base30cms040_all_Q base150cms040_all_Q

% Collect flow under steady peak flow for bank n = 0.040

peak180cms040_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3083_180cmsbank040/'...
    'abq10a12bf180cmsbank040rst02_0000_20160104_160000_his.nc'];
peak300cms040_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3084_300cmsbank040/'...
    'abq10a12bf300cmsbank040rst02_0000_20160103_210000_his.nc'];

% Extract flow at all cross-section times
peak180cms040_all_Q = ncread(peak180cms040_file,'cross_section_discharge');
peak300cms040_all_Q = ncread(peak300cms040_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms040_Q = peak180cms040_all_Q(1:max_dist/250*3,end);
peak300cms040_Q = peak300cms040_all_Q(1:max_dist/250*3,end);

clear peak180cms040_all_Q peak300cms040_all_Q

%% Collect peak flow information from bank n = 0.040 flood runs

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

% Define folder with *his.nc files
fpath = ['run3071_s180bank040/';'run3072_s300bank040/';...
    'run3073_m180bank040/';'run3074_m300bank040/';...
    'run3075_l180bank040/';'run3076_l300bank040/']; % file path

% fpath = ['runR025_bank040s180/';'runR026_bank040s300/';...
%     'runR027_bank040m180/';'runR028_bank040m300/';...
%     'runR029_bank040l180/';'runR030_bank040l300/']; % file path

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

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

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms040_Q;
        xs_base_Q = base150cms040_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms040_Q;
        xs_base_Q = base30cms040_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_bank040_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 n = 0.060
% Collect baseflow maximum for bank 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 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 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_bank060_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_bank060_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_bank060_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 n = 0.080
% Collect baseflow maximum for bank n = 0.080

base30cms080_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3069_30cmsbank100/'...
    'abq10a12bf30cmsbank100rst01_0000_20160104_090000_his.nc'];
base150cms080_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3070_150cmsbank100/'...
    'abq10a12bf150cmsbank100rst02_0000_20160104_210000_his.nc'];

% Extract flow at all cross-section times
base30cms080_all_Q = ncread(base30cms080_file,'cross_section_discharge');
base150cms080_all_Q = ncread(base150cms080_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms080_Q = base30cms080_all_Q(1:max_dist/250*3,end);
base150cms080_Q = base150cms080_all_Q(1:max_dist/250*3,end);

clear base30cms080_all_Q base150cms080_all_Q

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

peak180cms080_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3087_180cmsbank100/'...
    'abq10a12bf180cmsbank100rst02_0000_20160104_130000_his.nc'];
peak300cms080_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3088_300cmsbank100/'...
    'abq10a12bf300cmsbank100rst02_0000_20160103_200000_his.nc'];

% Extract flow at all cross-section times
peak180cms080_all_Q = ncread(peak180cms080_file,'cross_section_discharge');
peak300cms080_all_Q = ncread(peak300cms080_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms080_Q = peak180cms080_all_Q(1:max_dist/250*3,end);
peak300cms080_Q = peak300cms080_all_Q(1:max_dist/250*3,end);

clear peak180cms080_all_Q peak300cms080_all_Q

%% Collect peak flow information from bank n = 0.080 flood runs

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

% Define folder with *his.nc files
fpath = ['run3077_s180bank080/';'run3078_s300bank080/';...
    'run3079_m180bank080/';'run3080_m300bank080/';...
    'run3081_l180bank080/';'run3082_l300bank080/']; % file path

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

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

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms080_Q;
        xs_base_Q = base150cms080_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms080_Q;
        xs_base_Q = base30cms080_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_bank080_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 n = 0.100
% Collect baseflow maximum for bank n = 0.100

base30cms100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3067_30cmsbank080/'...
    'abq10a12bf30cmsbank080rst01_0000_20160104_080000_his.nc'];
base150cms100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3068_150cmsbank080/'...
    'abq10a12bf150cmsbank080rst01_0000_20160103_000000_his.nc'];

% Extract flow at all cross-section times
base30cms100_all_Q = ncread(base30cms100_file,'cross_section_discharge');
base150cms100_all_Q = ncread(base150cms100_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
base30cms100_Q = base30cms100_all_Q(1:max_dist/250*3,end);
base150cms100_Q = base150cms100_all_Q(1:max_dist/250*3,end);

clear base30cms100_all_Q base150cms100_all_Q

%% Collect flow under steady peak flow for bank n = 0.100

peak180cms100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3085_180cmsbank080/'...
    'abq10a12bf180cmsbank080rst02_0000_20160104_140000_his.nc'];
peak300cms100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/run3086_300cmsbank080/'...
    'abq10a12bf300cmsbank080rst02_0000_20160103_210000_his.nc'];

% Extract flow at all cross-section times
peak180cms100_all_Q = ncread(peak180cms100_file,'cross_section_discharge');
peak300cms100_all_Q = ncread(peak300cms100_file,'cross_section_discharge');

% Extract steady state (end time) flow at all cross-sections
peak180cms100_Q = peak180cms100_all_Q(1:max_dist/250*3,end);
peak300cms100_Q = peak300cms100_all_Q(1:max_dist/250*3,end);

clear peak180cms080_all_Q peak300cms080_all_Q

%% Collect peak flow information from bank n = 0.100 flood runs

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

% Define folder with *his.nc files
fpath = ['run3089_s180bank100/';'run3090_s300bank100/';...
    'run3091_m180bank100/';'run3092_m300bank100/';...
    'run3093_l180bank100/';'run3094_l300bank100/']; % file path

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

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

    if mod(fi,2) == 0
        xs_peak_Q = peak300cms100_Q;
        xs_base_Q = base150cms100_Q;
    elseif mod(fi,2) == 1
        xs_peak_Q = peak180cms100_Q;
        xs_base_Q = base30cms100_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_bank100_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
    
    fp040_attenuation = xs_fp040_cell{i,9};    
    fp060_attenuation = xs_fp060_cell{i,9};
    fp080_attenuation = xs_fp080_cell{i,9};
    fp100_attenuation = xs_fp100_cell{i,9};
    bank040_attenuation = xs_bank040_cell{i,9};    
    bank060_attenuation = xs_bank060_cell{i,9};
    bank080_attenuation = xs_bank080_cell{i,9};
    bank100_attenuation = xs_bank100_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),...
        fp040_attenuation(1:max_dist_id,1),'--');
    p3=plot(xs_dist(1:max_dist_id),...
        fp080_attenuation(1:max_dist_id,1),'--');
    p4=plot(xs_dist(1:max_dist_id),...
        fp100_attenuation(1:max_dist_id,1),'--');
    ax = gca;
    ax.ColorOrderIndex = 1;
    p5=plot(xs_dist,bank040_attenuation(1:max_dist/250,1));
    p6=plot(xs_dist,bank080_attenuation(1:max_dist/250,1));
    p7=plot(xs_dist,bank100_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 p4],'Baseline','F04','F08',...
%             'F10','Location','southwest');
%         legend('boxoff');
    elseif i == 6
        ylim([0.7 1]);
        xlabel('Downstream distance (km)');
%         legend([p5 p6 p7],'B04','B08','B10',...
%             'Location','southwest');
%         legend('boxoff');
    end
    xlim([0 28]);
    set(gca,'Fontsize',14);
    
    B04_end(i,1) = bank040_attenuation(max_dist/250,1);
    B08_end(i,1) = bank080_attenuation(max_dist/250,1);
    B10_end(i,1) = bank100_attenuation(max_dist/250,1);
    F04_end(i,1) = fp040_attenuation(max_dist/250,1);
    F08_end(i,1) = fp080_attenuation(max_dist/250,1);
    F10_end(i,1) = fp100_attenuation(max_dist/250,1);
    Baseline_end(i,1) = fp060_attenuation(max_dist/250,1);

    percent_B04(i,1) = (B04_end(i,1) - Baseline_end(i,1)) / ...
        (1 - Baseline_end(i,1));
    percent_B08(i,1) = (B08_end(i,1) - Baseline_end(i,1)) / ...
        (1 - Baseline_end(i,1));
    percent_B10(i,1) = (B10_end(i,1) - Baseline_end(i,1)) / ...
        (1 - Baseline_end(i,1));
    percent_F04(i,1) = (F04_end(i,1) - Baseline_end(i,1)) / ...
        (1 - Baseline_end(i,1));
    percent_F08(i,1) = (F08_end(i,1) - Baseline_end(i,1)) / ...
        (1 - Baseline_end(i,1));
    percent_F10(i,1) = (F10_end(i,1) - Baseline_end(i,1)) / ...
        (1 - Baseline_end(i,1));
    
end

max_topo_bank = max([max(percent_B04) max(percent_B08) max(percent_B10)]);
max_topo_ch = max([max(percent_F04) max(percent_F08) max(percent_F10)]);
min_topo_bank = min([min(percent_B04) min(percent_B08) min(percent_B10)]);
min_topo_ch = min([min(percent_F04) min(percent_F08) min(percent_F10)]);

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