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

%Indexing
J1=J(2:jamx-1);         K1=K(2:jmax-1);
J2=J(1:jmax-2);         K2=K(1:jmax-2);

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