% Prepares data from 1942 model runs and plots all in one plot

clear

%% Run *his analysis for each 1942 baseflow run

% Peak flow files

f1942_fnl_2010 = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012_histextent/run7023_f583cms/'...
    'abq10a12bshcf583cmsrst02_0000_20160104_060000_his.nc'];
f1942_fnl_hist = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abqhistorical/run6024_f583redo/'...
    'historicf583redorst03_0000_19180104_070000_his.nc'];
f1942_fnl_2010lidar = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12blhc/run8011_f583cms/'...
    'abq10a12blhcf583cmsrst02_0000_20160103_110000_his.nc'];
f1942_fnl_prelidar = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abqpreblhc/run9011_f583cms/'...
    'abqprebhcf583cmsrst02_0000_20160104_100000_his.nc'];

all_2010_f1942 = ncread(f1942_fnl_2010,'waterlevel');
all_hist_f1942 = ncread(f1942_fnl_hist,'waterlevel');
all_2010lidar_f1942 = ncread(f1942_fnl_2010lidar,'waterlevel');
all_prelidar_f1942 = ncread(f1942_fnl_prelidar,'waterlevel');

obs_2010_wse(:,1) = all_2010_f1942(:,end-1);
obs_hist_wse(:,1) = all_hist_f1942(:,end-1);
obs_2010lidar_wse(:,1) = all_2010lidar_f1942(:,end-1);
obs_prelidar_wse(:,1) = all_prelidar_f1942(:,end-1);

% Load local minimum elevations for contemporary and historical models
load(['/Users/Colin/Documents/MATLAB/DFlow/attenuation_chap1/data/'...
    'valley250minZcontemporary.mat']);
load(['/Users/Colin/Documents/MATLAB/DFlow/attenuation_chap1/data/'...
    'valley250minZhistorical.mat']);

%% Run *his analysis for each 1942 run

%% Historical _his.nc function call

chan_dist_max = 18000; % km
max_Q = 583.7; % m^3/s
min_t = 7200; % plots start at 2 hour
max_t = 96 * 3600; % hours * seconds in an hour - gives seconds
fpath = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abqhistorical/run6022_f1942full/'];

segment250pts = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'mass_Mo_chap3/abq_250m_segment_polygons'],' '); % import 250 segments

folder_num = numel(fpath(:,1));

xs_hist_cell = cell(folder_num,7);
obs_wse_hist_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_hist_cell(fi,1) = {fpath(fi,:)};
    obs_wse_hist_cell(fi,1) = {fpath(fi,:)};

    obs_peak_wse = obs_hist_wse(:,1);
    
    max_Q_fi = max_Q(fi,1);

    fnames = struct2cell(dir(strcat(fpath,'*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
    
    [xs_hist_cell(fi,2:7), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his(fnum, fnames, fpath, chan_dist_max, min_t, ...
        max_t, max_Q_fi);
    
    %% Extract observation point WSE information
    
    [obs_wse_hist_cell(fi,2:9)] = ...
        hpeak_obs_from_his(fnum, fnames, fpath, chan_dist_max, min_t,...
        max_t, obs_peak_wse,valley250minZhistorical);
    
end

% Area data retrieval

hist_583area = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'attenuation_chap1/data/historicf583redo_area.mat']);
hist_1942area = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'abq_models/abqhistorical/run6023_f1942area/'...
    'historicf1942area_area.mat']);

%% Contemporary contour _his.nc function call

chan_dist_max = 18000; % km
max_Q = 583.7; % m^3/s
min_t = 7200; % plots start at 2 hour
max_t = 96 * 3600; % hours * seconds in an hour - gives seconds
fpath = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012_histextent/run7015_1942daily/'];

segment250pts = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'mass_Mo_chap3/abq_250m_segment_polygons'],' '); % import 250 segments

folder_num = numel(fpath(:,1));

xs_2010_cell = cell(folder_num,7);
obs_wse_2010_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_2010_cell(fi,1) = {fpath(fi,:)};
    obs_wse_2010_cell(fi,1) = {fpath(fi,:)};

    obs_peak_wse = obs_2010_wse(:,1);
    
    max_Q_fi = max_Q(fi,1);
    
    fnames = struct2cell(dir(strcat(fpath,'*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
    
    [xs_2010_cell(fi,2:7), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his(fnum, fnames, fpath, chan_dist_max, min_t, ...
        max_t, max_Q_fi);
    
    %% Extract observation point WSE information
    
    [obs_wse_2010_cell(fi,2:9)] = ...
        hpeak_obs_from_his(fnum, fnames, fpath, chan_dist_max, min_t,...
        max_t, obs_peak_wse,valley250minZcontemporary);
    
end

% Area data retrieval

c2010_583area = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'abq_models/abq2010add2012_histextent/run7023_f583cms/'...
    'run7023_f583cms_area.mat']);
c2010_1942area = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'abq_models/abq2010add2012_histextent/run7016_1942dailyarea/'...
    'run7016_1942dailyarea_area.mat']);

%% Contemporary lidar _his.nc function call

chan_dist_max = 25000; % km
max_Q = 583.7; % m^3/s
min_t = 7200; % plots start at 2 hour
max_t = 96 * 3600; % hours * seconds in an hour - gives seconds
fpath = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq10a12blhc/run8019_f1942daily/'];

segment250pts = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'mass_Mo_chap3/abq_250m_segment_polygons'],' '); % import 250 segments

folder_num = numel(fpath(:,1));

xs_2010lidar_cell = cell(folder_num,7);
obs_wse_2010lidar_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_2010lidar_cell(fi,1) = {fpath(fi,:)};
    obs_wse_2010lidar_cell(fi,1) = {fpath(fi,:)};

    obs_peak_wse = obs_2010lidar_wse(:,1);
    
    max_Q_fi = max_Q(fi,1);
    
    fnames = struct2cell(dir(strcat(fpath,'*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
    
    [xs_2010lidar_cell(fi,2:7), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his(fnum, fnames, fpath, chan_dist_max, min_t, ...
        max_t, max_Q_fi);
    
    %% Extract observation point WSE information
    
    [obs_wse_2010lidar_cell(fi,2:9)] = ...
        hpeak_obs_from_his(fnum, fnames, fpath, chan_dist_max, min_t,...
        max_t, obs_peak_wse,valley250minZcontemporary);
    
end

% Area data retrieval

l2010_583area = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'abq_models/abq10a12blhc/run8011_f583cms/'...
    'run8011_f583cms_area.mat']);
l2010_1942area = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'abq_models/abq10a12blhc/run8020_f1942area/'...
    'run8020_f1942area_area.mat']);

%% Pre-restoration _his.nc function call

chan_dist_max = 25000; % km
max_Q = 583.7; % m^3/s
min_t = 7200; % plots start at 2 hour
max_t = 96 * 3600; % hours * seconds in an hour - gives seconds
fpath = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abqpreblhc/run9019_f1942daily/'];

segment250pts = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'mass_Mo_chap3/abq_250m_segment_polygons'],' '); % import 250 segments

folder_num = numel(fpath(:,1));

xs_prelidar_cell = cell(folder_num,7);
obs_wse_prelidar_cell = cell(folder_num,7);

for fi = 1:folder_num
    
    xs_prelidar_cell(fi,1) = {fpath(fi,:)};
    obs_wse_prelidar_cell(fi,1) = {fpath(fi,:)};

    obs_peak_wse = obs_prelidar_wse(:,1);
    
    max_Q_fi = max_Q(fi,1);
    
    fnames = struct2cell(dir(strcat(fpath,'*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
    
    [xs_prelidar_cell(fi,2:7), min_t_r, max_t_r, t_steps] = ...
        qpeak_xs_from_his(fnum, fnames, fpath, chan_dist_max, min_t, ...
        max_t, max_Q_fi);
    
    %% Extract observation point WSE information
    
    [obs_wse_prelidar_cell(fi,2:9)] = ...
        hpeak_obs_from_his(fnum, fnames, fpath, chan_dist_max, min_t,...
        max_t, obs_peak_wse,valley250minZcontemporary);
    
end

% Area data retrieval

p2010_583area = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'abq_models/abqpreblhc/run9011_f583cms/'...
    'run9011_f583cms_area.mat']);
p2010_1942area = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'abq_models/abqpreblhc/run9020_f1942area/'...
    'run9020_f1942area_area.mat']);

%% Creat figure

f1 = figure(1);
set(f1,'color','white');
fsize = 14;
h.Units = 'inches';
%h.Position = [1 1 11 15];
h.PaperSize = [8.5 11];

%% Plot 18-km hydrographs

xsid = 18*4;

sp1=subplot(4,1,1);
hist_time = xs_hist_cell{1,2};
hist_Q = xs_hist_cell{1,4};
hist_18kmQ = hist_Q(:,xsid);
c2010_time = xs_2010_cell{1,2};
c2010_Q = xs_2010_cell{1,4};
c2010_18kmQ = c2010_Q(:,xsid);
l2010_time = xs_2010lidar_cell{1,2};
l2010_Q = xs_2010lidar_cell{1,4};
l2010_18kmQ = l2010_Q(:,xsid);
p2010_time = xs_prelidar_cell{1,2};
p2010_Q = xs_prelidar_cell{1,4};
p2010_18kmQ = p2010_Q(:,xsid);

% Import upstream flood data
fld1942daily = importdata(['/Users/Colin/Documents/MATLAB/DFlow/'...
    'attenuation_chap1/data/fld1942daily.mat']);

plot(fld1942daily(:,1)/60/24,fld1942daily(:,2),'k-')
hold on
plot(hist_time/86400,hist_18kmQ)
plot(c2010_time/86400,c2010_18kmQ,'--')
plot(p2010_time/86400,p2010_18kmQ,'-.');
plot(l2010_time/86400,l2010_18kmQ,':');
xlim([0 6]);
xlabel('Time (days)');
ylabel({'Discharge'; '(m^3s^{-1})'});
legend('Upstream','HC','CC','PL','CL','Location','East');
set(gca,'Fontsize',fsize)

%% Plot QAR

sp2=subplot(4,1,2);
hist_dist = xs_hist_cell{1,6}/1000;
hist_QAR = xs_hist_cell{1,7};
c2010_dist = xs_2010_cell{1,6}/1000;
c2010_QAR = xs_2010_cell{1,7};
l2010_dist = xs_2010lidar_cell{1,6}/1000;
l2010_QAR = xs_2010lidar_cell{1,7};
p2010_dist = xs_prelidar_cell{1,6}/1000;
p2010_QAR = xs_prelidar_cell{1,7};

plot(hist_dist,hist_QAR)
hold on
plot(c2010_dist,c2010_QAR,'--')
plot(p2010_dist,p2010_QAR,'-.')
plot(l2010_dist,l2010_QAR,':')
xlabel('Distance downstream (km)')
ylabel('QAR')

set(gca,'Fontsize',fsize)

%% Plot SAR

sp3=subplot(4,1,3);
obs_dist_2010c = obs_wse_2010_cell{1,6}/1000;
obs_dist_hist = obs_wse_hist_cell{1,6}/1000;
obs_dist_2010l = obs_wse_2010lidar_cell{1,6}/1000;
obs_dist_2010p = obs_wse_prelidar_cell{1,6}/1000;

max_wse_hist = obs_wse_hist_cell{1,8};
max_wse_2010c = obs_wse_2010_cell{1,8};
max_wse_2010l = obs_wse_2010lidar_cell{1,8};
max_wse_2010p = obs_wse_prelidar_cell{1,8};


plot(obs_dist_hist,max_wse_hist)
hold on
plot(obs_dist_2010c,max_wse_2010c,'--')
plot(obs_dist_2010p,max_wse_2010p,'-.')
plot(obs_dist_2010l,max_wse_2010l,':')
xlabel('Distance downstream (km)')
ylabel('SAR')

set(gca,'Fontsize',fsize)

%% Plot IAR

sp4=subplot(4,1,4);

hist_iar = hist_1942area(1:(18*4),2) ./ hist_583area(1:(18*4),2);
c2010_iar = c2010_1942area(1:(18*4),2) ./ c2010_583area(1:(18*4),2);
l2010_iar = l2010_1942area(:,2) ./ l2010_583area(:,2);
p2010_iar = p2010_1942area(:,2) ./ p2010_583area(:,2);

plot(obs_dist_hist,hist_iar)
hold on
plot(obs_dist_2010c,c2010_iar,'--')
plot(obs_dist_2010p,p2010_iar,'-.')
plot(obs_dist_2010l,l2010_iar,':')
xlabel('Distance downstream (km)')
ylabel('IAR')

set(gca,'Fontsize',fsize)

%% Adjust locations of subplots

% sp1pos = get(sp1,'Position');
% sp2pos = get(sp2,'Position');
% sp3pos = get(sp3,'Position');
% sp4pos = get(sp4,'Position');
% 
% sp1pos(2) = 0.7673;
% sp1pos(4) = 0.1577;

%% Print to page

print('-fillpage',['/Users/Colin/Documents/MATLAB/DFlow/'...
    'attenuation_chap1/figures/f1942daily_v2_8_5x11'],'-dpdf');