function F = btrix_3x3(IL,IU,KM,B,C,D,R)
%
%            3X3 BLOCK Tri-DIAGONAL SOLVER
%
%USE HERE:  IL=2    IU=Jmax-1     KM=KM     B=b_x   C=c_x   D=d_x   R=R
IS = IL + 2;
I  = IL;
ierr = 0;
ior  = 0;
%
for n = 1:3
  for m = 1:3
    G(KM,n,m) = C(I,KM,n,m);
  end
end
[L11,L21,L22,L31,L32,L33,V1,V2,V3,U12,U13,U23,ierr] = ludecp(G,KM,ior,ierr);
if ~(ierr==0)
  fprintf('1:Failure of LU decomposition ierr = %g',ierr);
end

D1 = V1.*R(I,KM,1);
D2 = V2.*( R(I,KM,2) - L21.*D1);
D3 = V3.*( R(I,KM,3) - L31.*D1 - L32.*D2);
R(I,KM,3) = D3;
R(I,KM,2) = D2 - U23.*R(I,KM,3);
R(I,KM,1) = D1 - U13.*R(I,KM,3) - U12.*R(I,KM,2);
%
for M=1:3
  D1 = V1.*D(I,KM,1,M);
  D2 = V2.*(D(I,KM,2,M) - L21.*D1);
  D3 = V3.*(D(I,KM,3,M) - L31.*D1 - L32.*D2);
  U(I,KM,3,M) = D3;
  U(I,KM,2,M) = D2 - U23.*U(I,KM,3,M);
  U(I,KM,1,M) = D1 - U13.*U(I,KM,3,M) - U12.*U(I,KM,2,M);
end
%
%
I = IL + 1;
L = I  - 1;
%
%
for N=1:3
  for M=1:3
    G(KM,N,M) = C(I,KM,N,M)                         ...
	- B(I,KM,N,1).*U(L,KM,1,M)                  ...
	- B(I,KM,N,2).*U(L,KM,2,M)                  ...
	- B(I,KM,N,3).*U(L,KM,3,M);
  end
end
[L11,L21,L22,L31,L32,L33,V1,V2,V3,U12,U13,U23,ierr] = ludecp(G,KM,ior,ierr);
if ~(ierr==0)
  fprintf('2:Failure of LU decomposition ierr = %g',ierr);
end
for M=1:3
  R(I,KM,M) = R(I,KM,M)                        ...
      - B(I,KM,M,1).*R(L,KM,1)                 ...
      - B(I,KM,M,2).*R(L,KM,2)                 ...
      - B(I,KM,M,3).*R(L,KM,3);
end
D1 = V1.*R(I,KM,1);
D2 = V2.*( R(I,KM,2) - L21.*D1);
D3 = V3.*( R(I,KM,3) - L31.*D1 - L32.*D2);
R(I,KM,3) = D3;
R(I,KM,2) = D2 - U23.*R(I,KM,3);
R(I,KM,1) = D1 - U13.*R(I,KM,3) - U12.*R(I,KM,2);
for M=1:3
  %     
  D1 = V1.*D(I,KM,1,M);
  D2 = V2.*(D(I,KM,2,M) - L21.*D1);
  D3 = V3.*(D(I,KM,3,M) - L31.*D1 - L32.*D2);
  U(I,KM,3,M) = D3;
  U(I,KM,2,M) = D2 - U23.*U(I,KM,3,M);
  U(I,KM,1,M) = D1 - U13.*U(I,KM,3,M) - U12.*U(I,KM,2,M);
end
%
%  
for I=IS:IU
  L = I - 1;
  LL= I - 2;
  for N=1:3
    for M=1:3
      G(KM,N,M) = C(I,KM,N,M)                    ...
	  - B(I,KM,N,1).*U(L,KM,1,M)             ...
	  - B(I,KM,N,2).*U(L,KM,2,M)             ...
	  - B(I,KM,N,3).*U(L,KM,3,M);
    end
  end
  [L11,L21,L22,L31,L32,L33,V1,V2,V3,U12,U13,U23,ierr] = ludecp(G,KM,ior,ierr);
  if ~(ierr==0)
    fprintf('3:Failure of LU decomposition ierr = %g',ierr);
  end
  for M=1:3
    R(I,KM,M) = R(I,KM,M)                         ...
	- B(I,KM,M,1).*R(L,KM,1)                  ...
	- B(I,KM,M,2).*R(L,KM,2)                  ...
	- B(I,KM,M,3).*R(L,KM,3);                  
  end
  D1 = V1.*R(I,KM,1);
  D2 = V2.*( R(I,KM,2) - L21.*D1);
  D3 = V3.*( R(I,KM,3) - L31.*D1 - L32.*D2);
  R(I,KM,3) = D3;
  R(I,KM,2) = D2 - U23.*R(I,KM,3);
  R(I,KM,1) = D1 - U13.*R(I,KM,3) - U12.*R(I,KM,2);
  while ~(I == IU)
    for M=1:3
      D1 = V1.*D(I,KM,1,M);
      D2 = V2.*(D(I,KM,2,M) - L21.*D1);
      D3 = V3.*(D(I,KM,3,M) - L31.*D1 - L32.*D2);
      U(I,KM,3,M) = D3;
      U(I,KM,2,M) = D2 - U23.*U(I,KM,3,M);
      U(I,KM,1,M) = D1 - U13.*U(I,KM,3,M) - U12.*U(I,KM,2,M);
    end
    if I == IU-1, break, end
    I = IU;
  end
end
%
%    BACK SWEEP
%
I = IU - 1;
for M=1:3
  R(I,KM,M) = R(I,KM,M) - U(I,KM,M,1).*R(I+1,KM,1)         ...
     - U(I,KM,M,2).*R(I+1,KM,2) - U(I,KM,M,3).*R(I+1,KM,3);
end
for I = IU-2:-1:IL
  for M=1:3
    R(I,KM,M) = R(I,KM,M)                                 ...
	- U(I,KM,M,1).*R(I+1,KM,1)                        ...
	- U(I,KM,M,2).*R(I+1,KM,2)                        ...
	- U(I,KM,M,3).*R(I+1,KM,3);
  end
end
F = R;