function [obs_return] =  hpeak_obs_from_his_chap2(...
    file_num, filename, file_loc, chan_dist_max, min_t, max_t, ...
    obs_peak_wse, minimum_xs_z, obs_xs_ids, obs_rcl_ids)
% Function takes *_his.nc DFlow model output and extracts peak depth and
% water surface elevation

for i = 1:file_num
    fnameloc(i) = {strcat(file_loc, filename(i))};
end

% Load water level data
th_raw = cell(file_num,2);
td_raw = cell(file_num,2);

for i = 1:file_num    
    th_raw(i,1) = {ncread(char(fnameloc{i}), 'time')};
    th_raw(i,2) = {ncread(char(fnameloc{i}), 'waterlevel')};
    
    td_raw(i,1) = {ncread(char(fnameloc{i}), 'time')};
    td_raw(i,2) = {ncread(char(fnameloc{i}), 'Waterdepth')};
end
% 
for x = 1:file_num
    if x == 1
        % Create WSE matrix
        time_wse(:,1) = cell2mat(th_raw(x,1));
        wse = cell2mat(th_raw(x,2));
        [obs_count,~] = size(wse);
        time_wse(:,2:obs_count+1) = wse';
        th_hold = 0;
        [t_count, ~] = size(time_wse);
        % Create Depth matrix
        time_depth(:,1) = cell2mat(td_raw(x,1));
        depth = cell2mat(td_raw(x,2));
        [obs_count,~] = size(depth);
        time_depth(:,2:obs_count+1) = depth';
        td_hold = 0;
        [t_count_d, ~] = size(time_depth);
    elseif x > 1
        % Append WSE matrix
        clear th_hold
        th_hold(:,1) = cell2mat(th_raw(x,1));
        h_hold = cell2mat(th_raw(x,2));
        th_hold(:,2:obs_count+1) = h_hold';
        [tsteps_hold,~] = size(th_hold);
        for j = 1:tsteps_hold
            if th_hold(j,1) > time_wse(t_count,1)
                time_wse(t_count+1:t_count+tsteps_hold-j+1,:)...
                    = th_hold(j:tsteps_hold,:);
                [t_count, ~] = size(time_wse);
                break
            end
        end
        % Append Depth matrix
        clear td_hold
        td_hold(:,1) = cell2mat(td_raw(x,1));
        d_hold = cell2mat(td_raw(x,2));
        td_hold(:,2:obs_count+1) = d_hold';
        [tsteps_hold_d,~] = size(td_hold);
        for j = 1:tsteps_hold_d
            if td_hold(j,1) > time_depth(t_count_d,1)
                time_depth(t_count_d+1:t_count_d+tsteps_hold_d-j+1,:)...
                    = td_hold(j:tsteps_hold_d,:);
                [t_count_d, ~] = size(time_depth);
                break
            end
        end
            
    end
end

time_wse = abs(time_wse);
wse = time_wse(:,2:end);
time_depth = abs(time_depth);
depth = time_depth(:,2:end);
min_time = min(time_wse(:,1));
max_time = max(time_wse(:,1));

max_xs_id = max(obs_xs_ids(:,1));
wse_channel_avg = zeros(t_count,max_xs_id);
for i = 1:max_xs_id
    for j = 1:t_count
        wse_channel_avg(j,i) = mean(wse(j,obs_rcl_ids(:,1)'==2 &...
            obs_xs_ids(:,1)'==i & depth(j,:)>0));
    end
end

max_obs_num = chan_dist_max/250; % Extracts data from X km downstream
max_wse = zeros(max_obs_num, 3);
max_wse(:,1) = 1:max_obs_num;

for i = 1:max_obs_num

        max_wse(i,3) = max(wse_channel_avg(:,i));
        [~,r] = max(wse_channel_avg(:,i));
        max_wse(i,2) = time_wse(r,1);
        
end

obs_time = time_wse(:,1); % set time array
obs_wse_Xkm = wse_channel_avg(1:max_obs_num,:); % flow data for upper 27 km
obs_max_h_Xkm = max_wse(1:max_obs_num,:); % Extracts maximum h
obs_dist = obs_max_h_Xkm(:,1) * 250; % sets XS distances

% Compare to peak discharge model runs
obs_rarh = obs_max_h_Xkm(:,3) ./ obs_peak_wse(1:max_obs_num,1);
obs_rarh_minz = (obs_max_h_Xkm(:,3) - minimum_xs_z(1:max_obs_num,3)) ./...
    (obs_peak_wse(1:max_obs_num,1) - minimum_xs_z(1:max_obs_num,3));
% obs_rarh = (obs_peak_wse(1:max_obs_num,1) - obs_max_h_Xkm(:,3)) ./ ...
%     obs_peak_wse(1:max_obs_num,1); % calculates attenuation 

obs_return = cell(1,5);
obs_return(:,1) = {obs_time};
obs_return(:,2) = {obs_wse_Xkm};
obs_return(:,3) = {obs_max_h_Xkm};
obs_return(:,4) = {obs_dist};
obs_return(:,5) = {obs_rarh};
obs_return(:,6) = {obs_rarh_minz};

end