function [ RHSq ] = LS_NormalAdvect_UW(rho,NodeP,Nodeu,Nodev,dx,dy,idsu,...
                                       idsv,idsP,Nubcs,Nvbcs)
                                          
% LS_NORMALADVECT: 
% Advect Tracers (re-use RHSq to save a little memory)    
% To determine the RHS contribution of the level set expansion we solve
% phi_t + |grad phi| = 0 using upwind differences. So basically, we're
% looking for the upwind value of |grad phi|

[dpdx dpdy]   = Grad2D(rho,NodeP,dx,dy,idsu,idsv);
dpdx = [zeros(Nubcs,1);dpdx]; % Matrix of x derivatives of rho.
dpdy = [zeros(Nvbcs,1);dpdy]; % Matrix of y derivatives of rho.
%size(dpdx)
%size(dpdy)
%keyboard
tmpleft = [ones(size(NodeP, 1) - 2, 1),...
    rho(NodeP(2 : end - 1, 3 : end - 1))...
    <= rho(NodeP(2 : end - 1, 2 : end - 2))];
tmpright = [rho(NodeP(2 : end - 1, 2 : end - 2))...
    <= rho(NodeP(2 : end - 1, 3 : end - 1)),...
    ones(size(NodeP, 1)-2, 1)];
tmp =  ((1 - tmpleft) .* abs(...
    [dpdx(Nodeu(2 : end - 1, 2)), dpdx(Nodeu(2 : end - 1, 2 : end - 1))])...
    + (1- tmpright) .* abs(...
    [dpdx(Nodeu(2 : end - 1, 2 : end - 1)), dpdx(Nodeu(2 : end - 1, end - 1))]))...
    ./ max(tmpleft + tmpright, 1);
RHSq = tmp(idsP);

tmptop = [ones(1, size(NodeP, 2)-2);...
    rho(NodeP(3 : end - 1, 2 : end - 1))...
    <= rho(NodeP(2 : end - 2, 2 : end - 1))];
tmpbot = [rho(NodeP(2 : end - 2, 2 : end - 1))...
    <= rho(NodeP(3 : end - 1, 2 : end - 1));...
    ones(1, size(NodeP, 2)-2)];

tmp =  ((1 - tmptop) .* abs(...
    [dpdy(Nodev(2, 2 : end - 1))' ; dpdy(Nodev(2 : end - 1, 2 : end - 1))])...
    + (1 - tmpbot) .* abs(...
    [dpdy(Nodev(2 : end - 1, 2 : end - 1)) ; dpdy(Nodev(end - 1, 2 : end - 1))']))...
    ./ max(tmptop + tmpbot, 1);

RHSq = - sqrt(RHSq .^ 2 + (tmp(idsP)) .^ 2);

end

