function fc = AdvecSatur_FVM_QuadriQ0(EqnBIC,nx,ny,NumerFlux,sat)
%% Advection (convection only) of saturation in a 2-phase flow model 
% Solved by finite volumes on a logically rectangular quadrilateral mesh 
% James Liu, ColoState, Victor Ginting, UWYO; 2012/07--2018/12 

%% Fractional flow 
% fracfl = zeros(nx+2,ny+2);
fracfl = EqnBIC.fxnFractionalFlow(sat);

%% Fluxes on the left edges 
fluxLeft = zeros(nx+2,ny+2);
fluxLeft(2:nx+1,2:ny+1) = reshape(NumerFlux(:,4), nx, ny);
idx = fluxLeft<0.0;
ConcL = idx(2:nx+1,2:ny+1) .* fracfl(2:nx+1,1:ny)...
     + ~idx(2:nx+1,2:ny+1) .* fracfl(2:nx+1,2:ny+1);

%% Fluxes on the right edges 
fluxRight = zeros(nx+2,ny+2);
fluxRight(2:nx+1,2:ny+1) = reshape(NumerFlux(:,2), nx, ny);
idx = fluxRight<0.0;
ConcR = idx(2:nx+1,2:ny+1) .* fracfl(2:nx+1,3:ny+2)...
     + ~idx(2:nx+1,2:ny+1) .* fracfl(2:nx+1,2:ny+1);

%% Fluxes on the bottom edges 
fluxBottom = zeros(nx+2,ny+2);
fluxBottom(2:nx+1,2:ny+1) = reshape(NumerFlux(:,1), nx, ny);
idx = fluxBottom<0.0;
ConcB = idx(2:nx+1,2:ny+1) .* fracfl(1:nx,2:ny+1)...
     + ~idx(2:nx+1,2:ny+1) .* fracfl(2:nx+1,2:ny+1);

%% Fluxes on the top edges 
fluxTop = zeros(nx+2,ny+2);
fluxTop(2:nx+1,2:ny+1) = reshape(NumerFlux(:,3), nx, ny);
idx = fluxTop<0.0;
ConcT = idx(2:nx+1,2:ny+1) .* fracfl(3:nx+2,2:ny+1)...
     + ~idx(2:nx+1,2:ny+1) .* fracfl(2:nx+1,2:ny+1);

%% Summing up 
% fc = zeros(nx,ny);
fc = ConcL .* fluxLeft(2:nx+1,2:ny+1) ...
   + ConcR .* fluxRight(2:nx+1,2:ny+1) ...
   + ConcB .* fluxBottom(2:nx+1,2:ny+1) ...
   + ConcT .* fluxTop(2:nx+1,2:ny+1);

return;