function [obs_return] = ...
    hpeak_obs_from_his(file_num, filename, file_loc, chan_dist_max, ...
    min_t, max_t, obs_peak_wse, minimum_xs_z)
% 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 data from selected files for observation point location
obs_x = ncread(char(fnameloc{1}), 'station_x_coordinate');
obs_y = ncread(char(fnameloc{1}), 'station_y_coordinate');
obs_xy = [obs_x(:,1) obs_y(:,1)];

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
        time_wse(:,1) = cell2mat(th_raw(x,1));
        time_depth(:,1) = cell2mat(td_raw(x,1));
        wse = cell2mat(th_raw(x,2));
        depth = cell2mat(td_raw(x,2));
        [obs_count,~] = size(wse);
        time_wse(:,2:obs_count+1) = wse';
        time_depth(:,2:obs_count+1) = depth';
        th_hold = 0;
        td_hold = 0;
        [t_count, ~] = size(time_wse);
    elseif x > 1
        clear th_hold td_hold
        th_hold(:,1) = cell2mat(th_raw(x,1));
        td_hold(:,1) = cell2mat(td_raw(x,1));
        h_hold = cell2mat(th_raw(x,2));
        d_hold = cell2mat(td_raw(x,2));
        th_hold(:,2:obs_count+1) = h_hold';
        td_hold(:,2:obs_count+1) = d_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,:);
                time_depth(t_count+1:t_count+tsteps_hold-j+1,:)...
                    = td_hold(j:tsteps_hold,:);
                [t_count, ~] = size(time_wse);
                break
            end
        end
            
    end
end

time_wse = abs(time_wse);
time_depth = abs(time_depth);
min_time = min(time_wse(:,1));
max_time = max(time_wse(:,1));

max_wse = zeros(obs_count, 3);
max_wse(:,1) = 1:obs_count;

for i = 1:obs_count

        max_wse(i,3) = max(time_wse(:,i+1));
        [max_i,r] = max(time_wse(:,i+1));
        max_wse(i,2) = time_wse(r,1);
        %max_flows(i,4) = max_flows(i,3) / upstream_peak;
        
%         if i == obs_count
%             max_wse(i,1) = 17.5;
%         end
        
        %if i > 1
        %    max_flows(i,4) = 1000 / (max_flows(i,2) - max_flows(i-1,2));
        %else
        %    max_flows(i,4) = NaN;
        %end
        
end

max_obs_num = chan_dist_max/250; % Extracts data from X km downstream
obs_time = time_wse(:,1); % set time array
[min_t_r,~] = find(obs_time(:,1)==min_t); % index of set minimum time
[max_t_r,~] = find(obs_time(:,1)==max_t); % index of set maximum time
[t_steps,~] = size(obs_time); % total time steps in model run
obs_xy_Xkm = obs_xy(2:(max_obs_num+1),:); % Coordinates of XS
obs_flow_Xkm = time_wse(:, 2:(max_obs_num+1)); % h data for upper 27 km
obs_depth_Xkm = time_depth(:, 2:(max_obs_num+1)); % depth 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_xy_Xkm};
obs_return(:,3) = {obs_flow_Xkm};
obs_return(:,4) = {obs_max_h_Xkm};
obs_return(:,5) = {obs_dist};
obs_return(:,6) = {obs_rarh};
obs_return(:,7) = {obs_rarh_minz};
obs_return(:,8) = {obs_depth_Xkm};

end