function [R] = nl_rhs_cen(Q,JM,KM,dx,dy,epse)

%Get p,u, and v
p=Q(:,:,1);
u=Q(:,:,2)./Q(:,:,1);
v=Q(:,:,3)./Q(:,:,1);

sizeQ=size(Q);
N=sizeQ(1);
M=sizeQ(2);

A=zeros(N,M,3,3);
B=zeros(N,M,3,3);

A(:,:,1,1)=0;         A(:,:,1,2)=1;       A(:,:,1,3)=0;
A(:,:,2,1)=-u.^2+1;   A(:,:,2,2)=2*u;     A(:,:,2,3)=0;
A(:,:,3,1)=-u.*v;     A(:,:,3,2)=v;       A(:,:,3,3)=u;

B(:,:,1,1)=0;         B(:,:,1,2)=0;       B(:,:,1,3)=1;
B(:,:,2,1)=-u.*v;     B(:,:,2,2)=v;       B(:,:,2,3)=u;
B(:,:,3,1)=-v.^2;     B(:,:,3,2)=0;       B(:,:,3,3)=2*v;

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(JM,KM,m,n).*Qx_c(JM,KM,n) + ...
	B(JM,KM,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*Q(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,:);