function [R] = rhs_cen(Q,JM,KM,dx,dy,A,B,epse)

R = zeros(size(Q));
[Qx_c,Qy_c] = diff_cen(Q,JM,KM,dx,dy); 

for m = 1:3
  R(JM,KM,m) = 0.0;
  for n = 1:3
    R(JM,KM,m) = R(JM,KM,m) + ...
	A(m,n)*Qx_c(JM,KM,n) + ...
	B(m,n)*Qy_c(JM,KM,n) ;
  end
end

if epse ~= 0.0  
  Qx_d(JM,KM,:) = (Q(JM+1,KM,:) -2*Q(JM,KM,:) + Q(JM-1,KM,:))/dx;
  Qy_d(JM,KM,:) = (Q(JM,KM+1,:) -2*iQ(JM,KM,:) + Q(JM,KM-1,:))/dy;
  R(JM,KM,:) = R(JM,KM,:)-epse*(Qx_d(JM,KM,:) + Qy_d(JM,KM,:));
end

R(JM,KM,:) = -R(JM,KM,:);