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

% 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,:);
xs_elevation_up = xs_elevation;
xs_elevation_up(:,3) = xs_elevation_up(:,3) + 0.3048;
xs_elevation_down = xs_elevation;
xs_elevation_down(:,3) = xs_elevation_down(:,3) - 0.3048;

%% 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_wse = ncread(...
    peak180cms_chdown_file,'waterlevel');
peak300cms_chdown_all_wse = ncread(...
    peak300cms_chdown_file,'waterlevel');

% Extract steady state (end time) flow at all cross-sections
peak180cms_chdown_wse_steady = peak180cms_chdown_all_wse(:,end);
peak300cms_chdown_wse_steady = peak300cms_chdown_all_wse(:,end);

% Extract WSE at all cross-section times
peak180cms_chdown_all_depth = ncread(peak180cms_chdown_file,'Waterdepth');
peak300cms_chdown_all_depth = ncread(peak300cms_chdown_file,'Waterdepth');

% Extract steady state (end time) flow at all cross-sections
peak180cms_chdown_depth_steady = peak180cms_chdown_all_depth(:,end);
peak300cms_chdown_depth_steady = peak300cms_chdown_all_depth(:,end);

clear peak180cms_chdown_all_wse peak300cms_chdown_all_wse...
    peak180cms_chdown_all_depth peak300cms_chdown_all_depth

% Find maximum WSE for peak flow at each XS
peak180cms_chdown_wse = zeros(max_xs_num,1);
peak300cms_chdown_wse = zeros(max_xs_num,1);
for i = 1:max_xs_num
    peak180cms_chdown_wse(i,1) = mean(...
        peak180cms_chdown_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak180cms_chdown_depth_steady(:,1)>0,1));
    peak300cms_chdown_wse(i,1) = mean(...
        peak300cms_chdown_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak300cms_chdown_depth_steady(:,1)>0,1));
end

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

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

    if mod(fi,2) == 0
        obs_peak_wse = peak300cms_chdown_wse;
        %xs_base_Q = base150cms_chdown_Q;
    elseif mod(fi,2) == 1
        obs_peak_wse = peak180cms_chdown_wse;
        %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
    
    [obs_chdown_cell(fi,2:7)] = ...
        hpeak_obs_from_his_chap2(fnum, fnames, fpath_fi, max_dist,min_t,...
        max_t, obs_peak_wse, xs_elevation_down, obs_xs_ids, obs_rcl_ids);
    
end
disp('Done Channel Down')

%% 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_wse = ncread(...
    peak180cms_chup_file,'waterlevel');
peak300cms_chup_all_wse = ncread(...
    peak300cms_chup_file,'waterlevel');

% Extract steady state (end time) flow at all cross-sections
peak180cms_chup_wse_steady = peak180cms_chup_all_wse(:,end);
peak300cms_chup_wse_steady = peak300cms_chup_all_wse(:,end);

% Extract WSE at all cross-section times
peak180cms_chup_all_depth = ncread(peak180cms_chup_file,'Waterdepth');
peak300cms_chup_all_depth = ncread(peak300cms_chup_file,'Waterdepth');

% Extract steady state (end time) flow at all cross-sections
peak180cms_chup_depth_steady = peak180cms_chup_all_depth(:,end);
peak300cms_chup_depth_steady = peak300cms_chup_all_depth(:,end);

clear peak180cms_chup_all_wse peak300cms_chup_all_wse...
    peak180cms_chup_all_depth peak300cms_chup_all_depth

% Find maximum WSE for peak flow at each XS
peak180cms_chup_wse = zeros(max_xs_num,1);
peak300cms_chup_wse = zeros(max_xs_num,1);
for i = 1:max_xs_num
    peak180cms_chup_wse(i,1) = mean(...
        peak180cms_chup_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak180cms_chup_depth_steady(:,1)>0,1));
    peak300cms_chup_wse(i,1) = mean(...
        peak300cms_chup_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak300cms_chup_depth_steady(:,1)>0,1));
end

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

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

    if mod(fi,2) == 0
        obs_peak_wse = peak300cms_chup_wse;
        %xs_base_Q = base150cms_chup_Q;
    elseif mod(fi,2) == 1
        obs_peak_wse = peak180cms_chup_wse;
        %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
    
    [obs_chup_cell(fi,2:7)] = ...
        hpeak_obs_from_his_chap2(fnum, fnames, fpath_fi, max_dist,min_t,...
        max_t, obs_peak_wse, xs_elevation_up, obs_xs_ids, obs_rcl_ids);
    
end
disp('Done Channel Up')

%% 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_wse = ncread(...
%     base30cms_bankdown_file,'waterlevel');
% base150cms_bankdown_all_wse = ncread(...
%     base150cms_bankdown_file,'waterlevel');
% 
% % Extract steady state (end time) flow at all cross-sections
% base30cms_bankdown_wse = base30cms_bankdown_all_wse(:,end);
% base150cms_bankdown_wse = base150cms_bankdown_all_wse(:,end);

%% 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_wse = ncread(...
    peak180cms_bankdown_file,'waterlevel');
peak300cms_bankdown_all_wse = ncread(...
    peak300cms_bankdown_file,'waterlevel');

% Extract steady state (end time) flow at all cross-sections
peak180cms_bankdown_wse_steady = peak180cms_bankdown_all_wse(:,end);
peak300cms_bankdown_wse_steady = peak300cms_bankdown_all_wse(:,end);

% Extract WSE at all cross-section times
peak180cms_bankdown_all_depth = ...
    ncread(peak180cms_bankdown_file,'Waterdepth');
peak300cms_bankdown_all_depth = ...
    ncread(peak300cms_bankdown_file,'Waterdepth');

% Extract steady state (end time) flow at all cross-sections
peak180cms_bankdown_depth_steady = peak180cms_bankdown_all_depth(:,end);
peak300cms_bankdown_depth_steady = peak300cms_bankdown_all_depth(:,end);

clear base30cms_bankdown_all_wse base150cms_bankdown_all_wse...
    peak180cms_bankdown_all_depth peak300cms_bankdown_all_depth

% Find maximum WSE for peak flow at each XS
peak180cms_bankdown_wse = zeros(max_xs_num,1);
peak300cms_bankdown_wse = zeros(max_xs_num,1);
for i = 1:max_xs_num
    peak180cms_bankdown_wse(i,1) = mean(...
        peak180cms_bankdown_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak180cms_bankdown_depth_steady(:,1)>0,1));
    peak300cms_bankdown_wse(i,1) = mean(...
        peak300cms_bankdown_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak300cms_bankdown_depth_steady(:,1)>0,1));
end

%% 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
        obs_peak_wse = peak300cms_bankdown_wse;
        %xs_base_Q = base150cms_bankdown_wse;
    elseif mod(fi,2) == 1
        obs_peak_wse = peak180cms_bankdown_wse;
        %xs_base_Q = base30cms_bankdown_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_bankdown_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

%% 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(...
%     base150cms100_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_wse = ncread(...
    peak180cms_bankup_file,'waterlevel');
peak300cms_bankup_all_wse = ncread(...
    peak300cms_bankup_file,'waterlevel');

% Extract steady state (end time) flow at all cross-sections
peak180cms_bankup_wse_steady = peak180cms_bankup_all_wse(:,end);
peak300cms_bankup_wse_steady = peak300cms_bankup_all_wse(:,end);

% Extract WSE at all cross-section times
peak180cms_bankup_all_depth = ...
    ncread(peak180cms_bankup_file,'Waterdepth');
peak300cms_bankup_all_depth = ...
    ncread(peak300cms_bankup_file,'Waterdepth');

% Extract steady state (end time) flow at all cross-sections
peak180cms_bankup_depth_steady = peak180cms_bankup_all_depth(:,end);
peak300cms_bankup_depth_steady = peak300cms_bankup_all_depth(:,end);

clear peak180cms_bankup_all_wse peak300cms_bankup_all_wse...
    peak180cms_bankup_all_depth peak300cms_bankup_all_depth

% Find maximum WSE for peak flow at each XS
peak180cms_bankup_wse = zeros(max_xs_num,1);
peak300cms_bankup_wse = zeros(max_xs_num,1);
for i = 1:max_xs_num
    peak180cms_bankup_wse(i,1) = mean(...
        peak180cms_bankup_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak180cms_bankup_depth_steady(:,1)>0,1));
    peak300cms_bankup_wse(i,1) = mean(...
        peak300cms_bankup_wse_steady(obs_rcl_ids(:,1)==2 &...
        obs_xs_ids(:,1)==i & peak300cms_bankup_depth_steady(:,1)>0,1));
end

%% 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
        obs_peak_wse = peak300cms_bankup_wse;
        %xs_base_Q = base150cms_bankup_Q;
    elseif mod(fi,2) == 1
        obs_peak_wse = peak180cms_bankup_wse;
        %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
    
    [obs_bankup_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

%% 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_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/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));
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;
        %xs_base_Q = base150cms060_Q;
    elseif mod(fi,2) == 1
        obs_peak_wse = peak180cms060_wse;
        %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
    
    [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



%% Plot attenuation ratio
figure(2)
set(gcf,'color','white');
xs_dist = (250:250:max_dist)/1000;
for i = 1:6
        
    fp060_attenuation = obs_fp060_cell{i,6};
    chdown_attenuation = obs_chdown_cell{i,6};
    chup_attenuation = obs_chup_cell{i,6};
    bankdown_attenuation = xs_bankdown_cell{i,6};
    bankup_attenuation = xs_bankup_cell{i,6};
    
    max_dist_id = max_dist/250;
    
    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),chdown_attenuation(1:max_dist_id,1));
    plot(xs_dist(1:max_dist_id),chup_attenuation(1:max_dist_id,1));
    ax = gca;
    ax.ColorOrderIndex = 1;
    plot(xs_dist,bankdown_attenuation(1:max_dist/250,1),'--');
    plot(xs_dist,bankup_attenuation(1:max_dist/250,1),'--');
    if i == 1
        ylim([0.99985 1]);
    elseif i == 2
        ylim([0.99985 1]);
    elseif i == 3
        ylim([0.99985 1]);
        ylabel('Discharge Attenuation Ratio');
    elseif i == 4
        ylim([0.99985 1]);
    elseif i == 5 
        ylim([0.99985 1]);
        xlabel('Downstream distance (km)');
        legend('Baseline','CH Down','CH Up','Bank Down','Bank Up',...
            'Location','southwest');
    elseif i == 6
        ylim([0.99985 1]);
        xlabel('Downstream distance (km)');
    end
    xlim([0 28]);
    set(gca,'Fontsize',14);
    
end