%
% 2D Advection Diffusion Forward Model
%
% E. Di Lorenzo, March 2013
%


%========================================================
function RHS = ADVOFF_rhs(u,uvel,vvel,dx,dy,K,grd);
%========================================================

[I,J]=size(u);
rhs_diff = zeros(I,J);
rhs_adv = zeros(I,J);

% Load Velocity field
% NEED TO COMPLETE THIS TO UPLOAD PROPER TIME


up=uvel; up(up<0)=0;
um=uvel; um(um>0)=0;
vp=vvel; vp(vp<0)=0;
vm=vvel; vm(vm>0)=0;

x=2:I-1;
y=2:J-1;

% Diffusion Term
rhs_diff(x,y) = + K*( u(x+1,y) - 2*u(x,y) + u(x-1,y) )/dx^2 ...
    + K*( u(x,y+1) - 2*u(x,y) + u(x,y-1) )/dy^2  ;

% Select advection scheme
centered_scheme=0;  % 2nd order
upstream_scheme=1;  % 1st order
upstream_scheme2=0;  % 2nd order less diffusive


% Advection Term: upstream scheme
if upstream_scheme==1
    upls = u(x,y) -u(x-1,y  );      uneg = u(x+1,y  ) -u(x,y);
    vpls = u(x,y) -u(x  ,y-1);      vneg = u(x  ,y+1) -u(x,y);
    
    rhs_adv(x,y) = -(up(x,y).*upls + um(x,y).*uneg)/dx ...
        -(vp(x,y).*vpls + vm(x,y).*vneg)/dx;
end

% Advection Term: upstream scheme
if upstream_scheme2==1
    x=3:I-2;
    y=3:J-2;
    
    upls = 3*u(x,y) -4*u(x-1,y  ) +u(x-2,y  );      uneg = -u(x+2,y  ) + 4*u(x+1,y  ) -3*u(x,y);
    vpls = 3*u(x,y) -4*u(x  ,y-1) +u(x  ,y-2);      vneg = -u(x  ,y+2) + 4*u(x  ,y+1) -3*u(x,y);
    
    rhs_adv(x,y) = -(up(x,y).*upls + um(x,y).*uneg)/(2*dx) ...
        -(vp(x,y).*vpls + vm(x,y).*vneg)/(2*dy);
    
end


RHS=rhs_diff+rhs_adv;
return