function [F] = SCAadvect(Nx, Ny, dx, dy, P, u, v, NodeP, Nodeu, Nodev, idsP)
% function [F] = SCAadvect(Nx, Ny, dx, dy, P, u, v, NodeP, Nodeu, Nodev)
% This function calculates the FV fluxes in a staggered grid for Tracer
% advection using an upwind advection scheme
%
% INPUTS:
%       Nx:     Number of interior nodes in the x-direction
%       Ny:     Number of interior nodes in the y-direction
%       dx:     Step size in x direction
%       dy:     Step size in y direction
%       P:      Vector containing the values of Tracer, including boundary data
%       u:      Vector containing the u values of velocity, including boundary
%                 data
%       v:      Vector containing the v values of velocity, including boudary
%                 data
%       NodeP:  Node numbering matrix for the scalar
%       Nodeu:  Node numbering matrix for the u velocity
%       Nodev:  Node numbering matrix for the v velocity
%       idsP:   Location of active Tracer nodes in interior
% OUTPUTS:
%       F:      Output vector: (Fe-Fw)*dy/dV+(Fs-Fn)*dx/dV
%
% Authors:   Matt Ueckermann, Pierre Lermusiaux, Pat Haley for MIT course 2.29

%% Prelinimaries
%% Initialize variables
F  = zeros(Nx*Ny,1);
%Using upwind fluxes
%NOTE: There is an annoying thing where from n(orth) to s(south) index I 
%decreases, whereas in x-y coordinates, I would increase.

%% Calculate fluxes (Vectorized for speed)
%Fe
w=1;e=0;s=0;n=0;
F(:) = reshape(F,Ny,Nx) + 0.5*dy/(dx*dy)*(...
         u(Nodeu(2-n:end-1+s,2-w:end-1+e))  .* ( P(NodeP(2-n:end-1-n,2+e:end-1+e)) + P(NodeP(2+s:end-1+s,2-w:end-1-w)) )...
    -abs(u(Nodeu(2-n:end-1+s,2-w:end-1+e))) .* ( P(NodeP(2-n:end-1-n,2+e:end-1+e)) - P(NodeP(2+s:end-1+s,2-w:end-1-w)) )...
    );
%Fw
w=0;e=1;s=0;n=0;
F(:) = reshape(F,Ny,Nx) - 0.5*dy/(dx*dy)*(...
         u(Nodeu(2-n:end-1+s,2-w:end-1+e))  .* ( P(NodeP(2-n:end-1-n,2+e:end-1+e)) + P(NodeP(2+s:end-1+s,2-w:end-1-w)) )...
    -abs(u(Nodeu(2-n:end-1+s,2-w:end-1+e))) .* ( P(NodeP(2-n:end-1-n,2+e:end-1+e)) - P(NodeP(2+s:end-1+s,2-w:end-1-w)) )...
    );
%Fs
w=0;e=0;s=1;n=0;
F(:) = reshape(F,Ny,Nx) + 0.5*dx/(dx*dy)*(...
         v(Nodev(2-n:end-1+s,2-e:end-1+w))  .* ( P(NodeP(2-n:end-1-n,2+e:end-1+e)) + P(NodeP(2+s:end-1+s,2-w:end-1-w)) )...
    -abs(v(Nodev(2-n:end-1+s,2-e:end-1+w))) .* ( P(NodeP(2-n:end-1-n,2+e:end-1+e)) - P(NodeP(2+s:end-1+s,2-w:end-1-w)) )...
    );
%Fn
w=0;e=0;s=0;n=1;
F(:) = reshape(F,Ny,Nx) - 0.5*dx/(dx*dy)*(...
         v(Nodev(2-n:end-1+s,2-e:end-1+w))  .* ( P(NodeP(2-n:end-1-n,2+e:end-1+e)) + P(NodeP(2+s:end-1+s,2-w:end-1-w)) )...
    -abs(v(Nodev(2-n:end-1+s,2-e:end-1+w))) .* ( P(NodeP(2-n:end-1-n,2+e:end-1+e)) - P(NodeP(2+s:end-1+s,2-w:end-1-w)) )...
    );
F = F(idsP);
