% 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];
max_xs_num = max_dist/250+1;

% Load observation ids
id_floc = '/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/';
load(strcat(id_floc,'obs_xs_ids.mat'));
load(strcat(id_floc,'obs_rcl_ids.mat'));

% Load xs minimum elevations
load(strcat(id_floc,'data/rg2010_250m_min_z.mat'));
xs_elevation = rg2010_250m_min_z(2:end,:);

%% 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_wse = ncread(base30cms040_file,'waterlevel');
% base150cms040_all_wse = ncread(base150cms040_file,'waterlevel');
% 
% % Extract steady state (end time) flow at all cross-sections
% base30cms040_wse = base30cms040_all_wse(:,end);
% base150cms040_wse = base150cms040_all_wse(:,end);
% 
% clear base30cms040_all_wse base150cms040_all_wse

%% 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 WSE at all cross-section times
peak180cms040_all_wse = ncread(peak180cms040_file,'waterlevel');
peak300cms040_all_wse = ncread(peak300cms040_file,'waterlevel');

% Extract steady state (end time) flow at all cross-sections
peak180cms040_wse_steady = peak180cms040_all_wse(:,end);
peak300cms040_wse_steady = peak300cms040_all_wse(:,end);

% Extract WSE at all cross-section times
peak180cms040_all_depth = ncread(peak180cms040_file,'Waterdepth');
peak300cms040_all_depth = ncread(peak300cms040_file,'Waterdepth');

% Extract steady state (end time) flow at all cross-sections
peak180cms040_depth_steady = peak180cms040_all_depth(:,end);
peak300cms040_depth_steady = peak300cms040_all_depth(:,end);

clear peak180cms040_all_wse peak300cms040_all_wse ...
    peak180cms040_all_depth peak300cms040_all_depth

% Find maximum WSE for peak flow at each XS
peak180cms040_wse = zeros(max_xs_num,1);
peak300cms040_wse = zeros(max_xs_num,1);
for i = 1:max_xs_num
    peak180cms040_wse(i,1) = mean(...
        peak180cms040_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak180cms040_depth_steady(:,1)>0,1));
    peak300cms040_wse(i,1) = mean(...
        peak300cms040_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak300cms040_depth_steady(:,1)>0,1));
end

%% 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));
obs_fp040_cell = cell(folder_num,7);

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

    if mod(fi,2) == 0
        obs_peak_wse = peak300cms040_wse;
        %obs_base_wse = base150cms040_wse;
    elseif mod(fi,2) == 1
        obs_peak_wse = peak180cms040_wse;
        %obs_base_wse = base30cms040_wse;
    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
    
    [obs_fp040_cell(fi,2:7)] = ...
        hpeak_obs_from_his_chap2(fnum, fnames, fpath_fi, max_dist,min_t,...
        max_t, obs_peak_wse, xs_elevation, obs_xs_ids, obs_rcl_ids);
    
end
disp('Done FP 0.040')

%% 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/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(:,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_wse = ncread(...
    peak180cms060_file,'waterlevel');
peak300cms060_all_wse = ncread(...
    peak300cms060_file,'waterlevel');

% Extract steady state (end time) flow at all cross-sections
peak180cms060_wse_steady = peak180cms060_all_wse(:,end);
peak300cms060_wse_steady = peak300cms060_all_wse(:,end);

% Extract WSE at all cross-section times
peak180cms060_all_depth = ncread(peak180cms060_file,'Waterdepth');
peak300cms060_all_depth = ncread(peak300cms060_file,'Waterdepth');

% Extract steady state (end time) flow at all cross-sections
peak180cms060_depth_steady = peak180cms060_all_depth(:,end);
peak300cms060_depth_steady = peak300cms060_all_depth(:,end);

clear peak180cms060_all_wse peak300cms060_all_wse...
    peak180cms060_all_depth peak300cms060_all_depth

% Find maximum WSE for peak flow at each XS
peak180cms060_wse = zeros(max_xs_num,1);
peak300cms060_wse = zeros(max_xs_num,1);
for i = 1:max_xs_num
    peak180cms060_wse(i,1) = mean(...
        peak180cms060_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak180cms060_depth_steady(:,1)>0,1));
    peak300cms060_wse(i,1) = mean(...
        peak300cms060_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak300cms060_depth_steady(:,1)>0,1));
end

%% 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));
obs_fp060_cell = cell(folder_num,7);

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

    if mod(fi,2) == 0
        obs_peak_wse = peak300cms060_wse;
        %obs_base_wse = base150cms060_Q;
    elseif mod(fi,2) == 1
        obs_peak_wse = peak180cms060_wse;
        %obs_base_wse = 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
    
    [obs_fp060_cell(fi,2:7)] = ...
        hpeak_obs_from_his_chap2(fnum, fnames, fpath_fi, max_dist,min_t,...
        max_t, obs_peak_wse, xs_elevation, obs_xs_ids, obs_rcl_ids);
    
end
disp('Done FP 0.060')

%% 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_wse = ncread(...
    peak180cms080_file,'waterlevel');
peak300cms080_all_wse = ncread(...
    peak300cms080_file,'waterlevel');

% Extract steady state (end time) flow at all cross-sections
peak180cms080_wse_steady = peak180cms080_all_wse(:,end);
peak300cms080_wse_steady = peak300cms080_all_wse(:,end);

% Extract WSE at all cross-section times
peak180cms080_all_depth = ncread(peak180cms080_file,'Waterdepth');
peak300cms080_all_depth = ncread(peak300cms080_file,'Waterdepth');

% Extract steady state (end time) flow at all cross-sections
peak180cms080_depth_steady = peak180cms080_all_depth(:,end);
peak300cms080_depth_steady = peak300cms080_all_depth(:,end);

clear peak180cms080_all_wse peak300cms080_all_wse...
    peak180cms080_all_depth peak300cms080_all_depth

% Find maximum WSE for peak flow at each XS
peak180cms080_wse = zeros(max_xs_num,1);
peak300cms080_wse = zeros(max_xs_num,1);
for i = 1:max_xs_num
    peak180cms080_wse(i,1) = mean(...
        peak180cms080_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak180cms080_depth_steady(:,1)>0,1));
    peak300cms080_wse(i,1) = mean(...
        peak300cms080_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak300cms080_depth_steady(:,1)>0,1));
end

%% 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));
obs_fp080_cell = cell(folder_num,7);

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

    if mod(fi,2) == 0
        obs_peak_wse = peak300cms080_wse;
        %obs_base_wse = base150cms080_Q;
    elseif mod(fi,2) == 1
        obs_peak_wse = peak180cms080_wse;
        %obs_base_wse = 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
    
    [obs_fp080_cell(fi,2:7)] = ...
        hpeak_obs_from_his_chap2(fnum, fnames, fpath_fi, max_dist,min_t,...
        max_t, obs_peak_wse, xs_elevation, obs_xs_ids, obs_rcl_ids);
    
end
disp('Done FP 0.080')

%% 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_wse = ncread(...
    peak180cms100_file,'waterlevel');
peak300cms100_all_wse = ncread(...
    peak300cms100_file,'waterlevel');

% Extract steady state (end time) flow at all cross-sections
peak180cms100_wse_steady = peak180cms100_all_wse(:,end);
peak300cms100_wse_steady = peak300cms100_all_wse(:,end);

% Extract WSE at all cross-section times
peak180cms100_all_depth = ncread(peak180cms100_file,'Waterdepth');
peak300cms100_all_depth = ncread(peak300cms100_file,'Waterdepth');

% Extract steady state (end time) flow at all cross-sections
peak180cms100_depth_steady = peak180cms100_all_depth(:,end);
peak300cms100_depth_steady = peak300cms100_all_depth(:,end);

clear peak180cms100_all_wse peak300cms100_all_wse...
    peak180cms100_all_depth peak300cms100_all_depth

% Find maximum WSE for peak flow at each XS
peak180cms100_wse = zeros(max_xs_num,1);
peak300cms100_wse = zeros(max_xs_num,1);
for i = 1:max_xs_num
    peak180cms100_wse(i,1) = mean(...
        peak180cms100_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak180cms100_depth_steady(:,1)>0,1));
    peak300cms100_wse(i,1) = mean(...
        peak300cms100_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak300cms100_depth_steady(:,1)>0,1));
end

%% 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));
obs_fp100_cell = cell(folder_num,7);

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

    if mod(fi,2) == 0
        obs_peak_wse = peak300cms100_wse;
        %obs_base_wse = base150cms100_Q;
    elseif mod(fi,2) == 1
        obs_peak_wse = peak180cms100_wse;
        %obs_base_wse = 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
    
    [obs_fp100_cell(fi,2:7)] = ...
        hpeak_obs_from_his_chap2(fnum, fnames, fpath_fi, max_dist,min_t,...
        max_t, obs_peak_wse, xs_elevation, obs_xs_ids, obs_rcl_ids);
    
end
disp('Done FP 0.100')

%% 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_wse = ncread(base30cms040_file,'cross_section_discharge');
% base150cms040_all_wse = ncread(base150cms040_file,'cross_section_discharge');
% 
% % Extract steady state (end time) flow at all cross-sections
% base30cms040_wse = base30cms040_all_wse(1:max_dist/250*3,end);
% base150cms040_wse = base150cms040_all_wse(1:max_dist/250*3,end);
% 
% clear base30cms040_all_Q base150cms040_all_Q

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

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

% Extract flow at all cross-section times
peak180cmsb040_all_wse = ncread(...
    peak180cmsb040_file,'waterlevel');
peak300cmsb040_all_wse = ncread(...
    peak300cmsb040_file,'waterlevel');

% Extract steady state (end time) flow at all cross-sections
peak180cmsb040_wse_steady = peak180cmsb040_all_wse(:,end);
peak300cmsb040_wse_steady = peak300cmsb040_all_wse(:,end);

% Extract WSE at all cross-section times
peak180cmsb040_all_depth = ncread(peak180cmsb040_file,'Waterdepth');
peak300cmsb040_all_depth = ncread(peak300cmsb040_file,'Waterdepth');

% Extract steady state (end time) flow at all cross-sections
peak180cmsb040_depth_steady = peak180cmsb040_all_depth(:,end);
peak300cmsb040_depth_steady = peak300cmsb040_all_depth(:,end);

clear peak180cmsb040_all_wse peak300cmsb040_all_wse...
    peak180cmsb040_all_depth peak300cmsb040_all_depth

% Find maximum WSE for peak flow at each XS
peak180cmsb040_wse = zeros(max_xs_num,1);
peak300cmsb040_wse = zeros(max_xs_num,1);
for i = 1:max_xs_num
    peak180cmsb040_wse(i,1) = mean(...
        peak180cmsb040_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak180cmsb040_depth_steady(:,1)>0,1));
    peak300cmsb040_wse(i,1) = mean(...
        peak300cmsb040_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak300cmsb040_depth_steady(:,1)>0,1));
end

%% 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));
obs_bank040_cell = cell(folder_num,7);

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

    if mod(fi,2) == 0
        obs_peak_wse = peak300cmsb040_wse;
        %obs_base_wse = base150cms040_wse;
    elseif mod(fi,2) == 1
        obs_peak_wse = peak180cmsb040_wse;
        %obs_base_wse = base30cms040_wse;
    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
    
    [obs_bank040_cell(fi,2:7)] = ...
        hpeak_obs_from_his_chap2(fnum, fnames, fpath_fi, max_dist,min_t,...
        max_t, obs_peak_wse, xs_elevation, obs_xs_ids, obs_rcl_ids);
    
end
disp('Done Bank 0.040')

% %% 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/'...
%     'abq10a12bf180cmsfp060rst02_0000_20160104_220000_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_wse = ncread(peak180cms060_file,'cross_section_discharge');
% peak300cms060_all_wse = ncread(peak300cms060_file,'cross_section_discharge');
% 
% % Extract steady state (end time) flow at all cross-sections
% peak180cms060_wse_steady = peak180cms060_all_wse(1:max_dist/250*3,end);
% peak300cms060_wse_steady = peak300cms060_all_wse(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
%         obs_peak_wse = peak300cms060_wse_steady;
%         obs_base_wse = base150cms060_Q;
%     elseif mod(fi,2) == 1
%         obs_peak_wse = peak180cms060_wse_steady;
%         obs_base_wse = 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, obs_peak_wse, obs_base_wse, 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
% base30cmsb080_all_wse = ncread(base30cms080_file,'cross_section_discharge');
% base150cmsb080_all_wse = 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.080

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

% Extract flow at all cross-section times
peak180cmsb080_all_wse = ncread(...
    peak180cmsb080_file,'waterlevel');
peak300cmsb080_all_wse = ncread(...
    peak300cmsb080_file,'waterlevel');

% Extract steady state (end time) flow at all cross-sections
peak180cmsb080_wse_steady = peak180cmsb080_all_wse(:,end);
peak300cmsb080_wse_steady = peak300cmsb080_all_wse(:,end);

% Extract WSE at all cross-section times
peak180cmsb080_all_depth = ncread(peak180cmsb080_file,'Waterdepth');
peak300cmsb080_all_depth = ncread(peak300cmsb080_file,'Waterdepth');

% Extract steady state (end time) flow at all cross-sections
peak180cmsb080_depth_steady = peak180cmsb080_all_depth(:,end);
peak300cmsb080_depth_steady = peak300cmsb080_all_depth(:,end);

clear peak180cmsb080_all_wse peak300cmsb080_all_wse...
    peak180cmsb080_all_depth peak300cmsb080_all_depth

% Find maximum WSE for peak flow at each XS
peak180cmsb080_wse = zeros(max_xs_num,1);
peak300cmsb080_wse = zeros(max_xs_num,1);
for i = 1:max_xs_num
    peak180cmsb080_wse(i,1) = mean(...
        peak180cmsb080_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak180cmsb080_depth_steady(:,1)>0,1));
    peak300cmsb080_wse(i,1) = mean(...
        peak300cmsb080_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak300cmsb080_depth_steady(:,1)>0,1));
end

%% 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));
obs_bank080_cell = cell(folder_num,7);

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

    if mod(fi,2) == 0
        obs_peak_wse = peak300cmsb080_wse;
        %obs_base_wse = base150cms080_Q;
    elseif mod(fi,2) == 1
        obs_peak_wse = peak180cmsb080_wse;
        %obs_base_wse = 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
    
    [obs_bank080_cell(fi,2:7)] = ...
        hpeak_obs_from_his_chap2(fnum, fnames, fpath_fi, max_dist,min_t,...
        max_t, obs_peak_wse, xs_elevation, obs_xs_ids, obs_rcl_ids);
    
end
disp('Done Bank 0.080')

%% Bank n = 0.100
% Collect baseflow maximum for bank n = 0.100
% 
% base30cmsb100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
%     'abq2010add2012bath/run3067_30cmsbank080/'...
%     'abq10a12bf30cmsbank080rst01_0000_20160104_080000_his.nc'];
% base150cmsb100_file = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
%     'abq2010add2012bath/run3068_150cmsbank080/'...
%     'abq10a12bf150cmsbank080rst01_0000_20160103_000000_his.nc'];
% 
% % Extract flow at all cross-section times
% base30cmsb100_all_Q = ncread(base30cmsb100_file,'cross_section_discharge');
% base150cmsb100_all_Q = ncread(base150cmsb100_file,'cross_section_discharge');
% 
% % Extract steady state (end time) flow at all cross-sections
% base30cmsb100_Q = base30cmsb100_all_Q(1:max_dist/250*3,end);
% base150cmsb100_Q = base150cmsb100_all_Q(1:max_dist/250*3,end);
% 
% clear base30cmsb100_all_Q base150cmsb100_all_Q

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

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

% Extract flow at all cross-section times
peak180cmsb100_all_wse = ncread(...
    peak180cmsb100_file,'waterlevel');
peak300cmsb100_all_wse = ncread(...
    peak300cmsb100_file,'waterlevel');

% Extract steady state (end time) flow at all cross-sections
peak180cmsb100_wse_steady = peak180cmsb100_all_wse(:,end);
peak300cmsb100_wse_steady = peak300cmsb100_all_wse(:,end);

% Extract WSE at all cross-section times
peak180cmsb100_all_depth = ncread(peak180cmsb100_file,'Waterdepth');
peak300cmsb100_all_depth = ncread(peak300cmsb100_file,'Waterdepth');

% Extract steady state (end time) flow at all cross-sections
peak180cmsb100_depth_steady = peak180cmsb100_all_depth(:,end);
peak300cmsb100_depth_steady = peak300cmsb100_all_depth(:,end);

clear peak180cmsb100_all_wse peak300cmsb100_all_wse...
    peak180cmsb100_all_depth peak300cmsb100_all_depth

% Find maximum WSE for peak flow at each XS
peak180cmsb100_wse = zeros(max_xs_num,1);
peak300cmsb100_wse = zeros(max_xs_num,1);
for i = 1:max_xs_num
    peak180cmsb100_wse(i,1) = mean(...
        peak180cmsb100_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak180cmsb100_depth_steady(:,1)>0,1));
    peak300cmsb100_wse(i,1) = mean(...
        peak300cmsb100_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak300cmsb100_depth_steady(:,1)>0,1));
end

%% 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));
obs_bank100_cell = cell(folder_num,7);

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

    if mod(fi,2) == 0
        obs_peak_wse = peak300cmsb100_wse;
        %obs_base_wse = base150cms100_Q;
    elseif mod(fi,2) == 1
        obs_peak_wse = peak180cmsb100_wse;
        %obs_base_wse = 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
    
    [obs_bank100_cell(fi,2:7)] = ...
        hpeak_obs_from_his_chap2(fnum, fnames, fpath_fi, max_dist,min_t,...
        max_t, obs_peak_wse, xs_elevation, obs_xs_ids, obs_rcl_ids);
    
end
disp('Done Bank 0.100')

%% Plot attenuation ratio
figure(2)
set(gcf,'color','white');
xs_dist = (0:250:max_dist)/1000;
for i = 1:6
    
    fp040_attenuation = obs_fp040_cell{i,6};    
    fp060_attenuation = obs_fp060_cell{i,6};
    fp080_attenuation = obs_fp080_cell{i,6};
    fp100_attenuation = obs_fp100_cell{i,6};
    bank040_attenuation = obs_bank040_cell{i,6};    
    bank080_attenuation = obs_bank080_cell{i,6};
    bank100_attenuation = obs_bank100_cell{i,6};
    
    max_dist_id = max_dist/250 + 1;
    
    subplot(3,2,i);
    plot(xs_dist(1:max_dist_id),fp060_attenuation(1:max_dist_id,1),'k-');
    hold on
    plot(xs_dist(1:max_dist_id),fp040_attenuation(1:max_dist_id,1),'--');
    plot(xs_dist(1:max_dist_id),fp080_attenuation(1:max_dist_id,1),'--');
    plot(xs_dist(1:max_dist_id),fp100_attenuation(1:max_dist_id,1),'--');
    ax = gca;
    ax.ColorOrderIndex = 1;
    plot(xs_dist,bank040_attenuation(1:max_dist/250,1));
    plot(xs_dist,bank080_attenuation(1:max_dist/250,1));
    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.8 1]);
        ylabel('Stage Attenuation Ratio');
    elseif i == 4
        %ylim([0.8 1]);
    elseif i == 5 
        %ylim([0.9 1]);
        xlabel('Downstream distance (km)');
        legend('Baseline FP 0.060','FP 0.040','FP 0.080','FP 0.100',...
            'B 0.040','B 0.80','B 0.100','Location','southwest');
    elseif i == 6
        %ylim([0.9 1]);
        xlabel('Downstream distance (km)');
    end
    xlim([0 28]);
    set(gca,'Fontsize',14);
    
end