function [Q_transfer, Mo_flux_per_area, Mo_flux, Mo_unit_length,...
    area, length]= ...
    dflow_mass_momentum_his_calc(boundary_type,Vx_elems, Vy_elems,...
    depth_elems, time_steps, unit_vectors)
%Calculates mass and momentum transfer using DFlow data
%   Detailed explanation goes here

% Use nodes_xy to calculate vectors for Vx and Vy associated with nodes

rho = 1000; % kg/m^3
t = time_steps;
[r,~] = size(unit_vectors);
depth = zeros(r,t);
Vx = zeros(r,t);
Vy = zeros(r,t);
area = zeros(r,t);
length = zeros(r,1);
V_mag = zeros(r,t);
V_dot_n = zeros(r,t);
Q_transfer = zeros(r,t);
Mo_flux_per_area = zeros(r,t);
Mo_flux = zeros(r,t);
Mo_unit_length = zeros(r,t);

for m = 1:t
    for n = 1:r

        u = unit_vectors(n,3);
        v = unit_vectors(n,4);
        length(n,1) = unit_vectors(n,5);
        depth(n,m) = (depth_elems(n,m) + depth_elems(n+1,m)) / 2;
        Vx(n,m) = (Vx_elems(n,m) + Vx_elems(n+1,m)) / 2;
        Vy(n,m) = (Vy_elems(n,m) + Vy_elems(n+1,m)) / 2;
        area(n,m) = length(n,1) * depth(n,m); % m^2
        V_mag(n,m) = sqrt(Vx(n,m)^2 + Vy(n,m)^2); % m/s
        V_dot_n(n,m) = Vx(n,m) * u + Vy(n,m) * v; % m/s normal to floodplain boundary

        if (boundary_type == 0) || (boundary_type == 1)

            Q_transfer(n,m) = V_dot_n(n,m) * area(n,m); % m^3/s

            Mo_flux_per_area(n,m) = rho * V_dot_n(n,m) * abs(V_dot_n(n,m)); % N/m^2
            Mo_flux(n,m) = rho * V_dot_n(n,m) * abs(V_dot_n(n,m))...
                * area(n,m); % N
            Mo_unit_length(n,m) = Mo_flux(n,m) / length(n,1); % N / m

        elseif boundary_type == 2

            Q_transfer(n) = V_mag * area; % m^3/s

            Mo_flux_per_area(n) = rho * V_mag^2; % N/m^2
            Mo_flux(n) = rho * area * V_mag^2; % N

        end

    end
end

end

