function[R]=rhs_flux_split_so(Q,J,K,jmax,kmax,dx,dy,Ap,Am,Bp,Bm,epse)

JM=[2:jmax-1];         KM=[2:kmax-1];

R=zeros(size(Q));
[Qx_f,Qy_f]=so_diff_forward(Q,J,K,jmax,kmax,dx,dy);
[Qx_b,Qy_b]=so_diff_backward(Q,J,K,jmax,kmax,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,:);
 