function [Mo_return] = Mo_calcs(lb_cell, rb_cell, t_steps)
% Function organizes and calculates various momentum results along the
% channel-floodplain interface given raw *his.nc extraction files

% Load momentum flux
lb_Moflux = cell2mat(lb_cell(8,1));
rb_Moflux = cell2mat(rb_cell(8,1));

% Load length of interface elements
lb_length = cell2mat(lb_cell(13,1));
rb_length = cell2mat(rb_cell(13,1));

% Calculate number of elements on each interface
num_lb = numel(lb_length(:,1));
num_rb = numel(rb_length(:,1));

% Assign base momentum transfer
lb_Moflux_base = lb_Moflux(:,2);
rb_Moflux_base = rb_Moflux(:,2);

% Calculate rate of momentum transfer crossing boundary at each time step
lb_Moflux_no_base = zeros(num_lb,t_steps);
rb_Moflux_no_base = zeros(num_rb,t_steps);
lb_Moflux_no_base_per_l = zeros(num_lb,t_steps);
rb_Moflux_no_base_per_l = zeros(num_rb,t_steps);

for mo_i = 1:num_lb
    for mo_j = 1:t_steps
        lb_Moflux_no_base(mo_i,mo_j) = lb_Moflux(mo_i,mo_j) -...
            lb_Moflux_base(mo_i,1);
        lb_Moflux_no_base_per_l(mo_i,mo_j) = ...
            lb_Moflux_no_base(mo_i,mo_j) / lb_length(mo_i,1);
    end
end

for mo_i = 1:num_rb
    for mo_j = 1:t_steps
        rb_Moflux_no_base(mo_i,mo_j) = rb_Moflux(mo_i,mo_j) -...
            rb_Moflux_base(mo_i,1);
        rb_Moflux_no_base_per_l(mo_i,mo_j) =...
            rb_Moflux_no_base(mo_i,mo_j) / rb_length(mo_i,1);
    end
end

lb_Moflux_no_base(:,1) = 0;
rb_Moflux_no_base(:,1) = 0;
lb_Moflux_no_base_per_l(:,1) = 0;
rb_Moflux_no_base_per_l(:,1) = 0;

lb_Moflux_t = sum(lb_Moflux_no_base);
rb_Moflux_t = sum(rb_Moflux_no_base);

% Calculate positive and negative flux at left bank fp interface for
% each time step
lb_Mo_pos = zeros(num_lb,t_steps);
lb_Mo_neg = zeros(num_lb,t_steps);

for Mo_ii = 1:num_lb
    for Mo_jj = 1:t_steps
        if lb_Moflux_no_base(Mo_ii,Mo_jj) >= 0
            lb_Mo_pos(Mo_ii,Mo_jj) = lb_Moflux_no_base(Mo_ii,Mo_jj);
            lb_Mo_neg(Mo_ii,Mo_jj) = 0;
        elseif lb_Moflux_no_base(Mo_ii,Mo_jj) < 0
            lb_Mo_pos(Mo_ii,Mo_jj) = 0;
            lb_Mo_neg(Mo_ii,Mo_jj) = lb_Moflux_no_base(Mo_ii,Mo_jj);
        end
    end
end

lb_Mo_t_pos = sum(lb_Mo_pos);
lb_Mo_t_neg = sum(lb_Mo_neg);

% Calculate positive and negative flux at right bank fp interface for each
% time step
rb_Mo_pos = zeros(num_rb,t_steps);
rb_Mo_neg = zeros(num_rb,t_steps);

for Mo_ii = 1:num_rb
    for Mo_jj = 1:t_steps
        if rb_Moflux_no_base(Mo_ii,Mo_jj) >= 0
            rb_Mo_pos(Mo_ii,Mo_jj) = rb_Moflux_no_base(Mo_ii,Mo_jj);
            rb_Mo_neg(Mo_ii,Mo_jj) = 0;
        elseif rb_Moflux_no_base(Mo_ii,Mo_jj) < 0
            rb_Mo_pos(Mo_ii,Mo_jj) = 0;
            rb_Mo_neg(Mo_ii,Mo_jj) = rb_Moflux_no_base(Mo_ii,Mo_jj);
        end
    end
end

rb_Mo_t_pos = sum(rb_Mo_pos);
rb_Mo_t_neg = sum(rb_Mo_neg);

Mo_return = cell(1,12);
Mo_return(1,1) = {lb_Moflux};
Mo_return(1,2) = {rb_Moflux};
Mo_return(1,3) = {lb_Moflux_base};
Mo_return(1,4) = {rb_Moflux_base};
Mo_return(1,5) = {lb_Moflux_no_base};
Mo_return(1,6) = {rb_Moflux_no_base};
Mo_return(1,7) = {lb_Moflux_no_base_per_l};
Mo_return(1,8) = {rb_Moflux_no_base_per_l};
Mo_return(1,9) = {lb_Moflux_t};
Mo_return(1,10) = {rb_Moflux_t};
Mo_return(1,11) = {lb_Mo_t_pos};
Mo_return(1,12) = {lb_Mo_t_neg};
Mo_return(1,13) = {rb_Mo_t_pos};
Mo_return(1,14) = {rb_Mo_t_neg};

end