function [F NBids] = Narrowband_LS_TVD_mex(Nx, Ny, dx, dy, dt, rho, u, v, NodeP, Nodeu, Nodev, idsP,Nbcs,NbcsDP,NbcsDu,NbcsDv,NBlen)
% 
% 
% This function calculates the FV fluxes in a staggered grid for Level set
% advection using TVD advection scheme in a narrow band around the zero
%level set.
%
% 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
%       dt:     Time step size
%       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
%       Nbcs:   Number of boundary conditions
%       NbcsDP: Number of Dirichlet Tracer boundary conditions
%       NbcsDu: Number of Dirichlet u-velocity boundary conditions
%       NbcsDv: Number of Dirichlet v-velocity boundary conditions
% OUTPUTS:
%       F:      Output vector: (Fe-Fw)*dy/dV+(Fs-Fn)*dx/dV
%
% Authors:   Matt Ueckermann, Pierre Lermusiaux, Pat Haley for MIT course 2.29

%% Call C++ function to do advection. Note the casting to int32's (important)
[F NBids] = Narrowband_LS_mx(int32(Nx),int32(Ny),dx,dy,dt,rho,u,v,...
    int32(NodeP),int32(Nodeu),int32(Nodev),...
    int32(Nbcs),int32(NbcsDP),int32(NbcsDu),int32(NbcsDv),NBlen);
NBids = logical(NBids);
F = F(idsP);
