function [Q_return] = Q_calcs(lb_cell, rb_cell, t_steps)
% Function organizes and calculates discharge extracted along
% channel-floodplain interface

% Extract raw discharge at each element
lb_Q = cell2mat(lb_cell(6,1));
rb_Q = cell2mat(rb_cell(6,1));

% Get size of left and right boundaries
[num_lb,~] = size(lb_Q);
[num_rb,~] = size(rb_Q);

% Adjust starting flood wave discharge to baseflow discharge
lb_Q_above_base = zeros(num_lb,t_steps);
rb_Q_above_base = zeros(num_rb,t_steps);

% Set baseflow to the second time step as the first time step can be zero
lb_base_Q = lb_Q(:,12);
rb_base_Q = rb_Q(:,12);

for lb_i = 1:num_lb
    for lb_t = 1:t_steps
        % Uses restart value (i.e. the value at time step = 2 because time
        % step 1 is equal to zero everywhere), for flood waves, restart
        % value begins at a base flow rst file
        lb_Q_above_base(lb_i,lb_t) = lb_Q(lb_i,lb_t) - lb_base_Q(lb_i,1);
        if lb_Q_above_base(lb_i,lb_t) < 0.01 &&...
                lb_Q_above_base(lb_i,lb_t) > -0.01
            lb_Q_above_base(lb_i,lb_t) = 0;
        end
    end
end

lb_Q_above_base(:,1) = 0;

for rb_i = 1:num_rb
    for rb_t = 1:t_steps
        % Uses restart value (i.e. the value at time step = 2 because time
        % step 1 is equal to zero everywhere), for flood waves, restart
        % value begins at a base flow rst file        
        rb_Q_above_base(rb_i,rb_t) = rb_Q(rb_i,rb_t) - rb_base_Q(rb_i,1);
        if rb_Q_above_base(rb_i,rb_t) < 0.01 &&...
                rb_Q_above_base(rb_i,rb_t) > -0.01
            rb_Q_above_base(rb_i,rb_t) = 0;
        end
    end
end

rb_Q_above_base(:,1) = 0;

% Calculate discharge crossing boundary at each time step
lb_Q_t = sum(lb_Q_above_base);
rb_Q_t = sum(rb_Q_above_base);

Q_return = cell(1,6);
Q_return(1,1) = {lb_Q};
Q_return(1,2) = {rb_Q};
Q_return(1,3) = {lb_base_Q};
Q_return(1,4) = {rb_base_Q};
Q_return(1,5) = {lb_Q_above_base};
Q_return(1,6) = {rb_Q_above_base};
Q_return(1,7) = {lb_Q_t};
Q_return(1,8) = {rb_Q_t};

end