function[Qx_f,Qy_f]=diff_forward(Q,JM,KM,dx,dy)

Qx_f(JM,KM,:)=(Q(JM+1,KM,:)-Q(JM,KM,:))/dx;
Qy_f(JM,KM,:)=(Q(JM,KM+1,:)-Q(JM,KM,:))/dy;

