function[Qx_f,Qy_f]=so_diff_forward(Q,J,K,jmax,kmax,dx,dy)

%Indexing
J1=[2:jmax-1];         K1=[2:kmax-1];
J2=[1:jmax-2];         K2=[1:kmax-2];

%Qx_f(J1,K1,1:3)=0;
%Qy_f(J1,K1,1:3)=0;
Qx_f=zeros(size(Q));
Qy_f=zeros(size(Q));

%Difference Interior Points:                 2nd order scheme
Qx_f(J2,K1,:)=.5*(-3*Q(J2,K1,:)+4*Q(J2+1,K1,:)-Q(J2+2,K1,:))/dx;
Qy_f(J1,K2,:)=.5*(-3*Q(J1,K2,:)+4*Q(J1,K2+1,:)-Q(J1,K2+2,:))/dy;

%Difference Right and Top boundaries:        1st order scheme
Qx_f(jmax-1,K1,:)=(Q(jmax,K1,:)-Q(jmax-1,K1,:))/dx;
Qx_f(J1,kmax-1,:)=(Q(J1,kmax,:)-Q(J1,kmax-1,:))/dy;