cd function[R]=rhs_flux_split(Q,JM,KM,dx,dy,Ap,Am,Bp,Bm,epse)

R=zeros(size(Q));
[Qx_f,Qy_f]=diff_forward(Q,JM,KM,dx,dy);
[Qx_b,Qy_b]=diff_backward(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)...
            -Ap(m,n)*Qx_b(JM,KM,n)-Am(m,n)*Qx_f(JM,KM,n)...
            -Bp(m,n)*Qy_b(JM,KM,n)-Bm(m,n)*Qy_f(JM,KM,n);
     end
 end
 
 R(JM,KM,:)=-R(JM,KM,:);