anchor.diff:        Dlu0(9*i  )=Dlu0(9*i  )*SIGMA_DIAG
       enddo
   
  +    i = 1
  +    call ILU1a33 (DkINV, &
  +         Dlu0(9*i-8), Dlu0(9*i-7), Dlu0(9*i-6), &
  +         Dlu0(9*i-5), Dlu0(9*i-4), Dlu0(9*i-3), &
  +         Dlu0(9*i-2), Dlu0(9*i-1), Dlu0(9*i  ))
  +    Dlu0(9*i-8)= DkINV(1,1)
  +    Dlu0(9*i-7)= DkINV(1,2)
  +    Dlu0(9*i-6)= DkINV(1,3)
  +    Dlu0(9*i-5)= DkINV(2,1)
  +    Dlu0(9*i-4)= DkINV(2,2)
  +    Dlu0(9*i-3)= DkINV(2,3)
  +    Dlu0(9*i-2)= DkINV(3,1)
  +    Dlu0(9*i-1)= DkINV(3,2)
  +    Dlu0(9*i  )= DkINV(3,3)
  +
       do i= 2, NP
         IW1= 0
         IW2= 0
  
positive.code:     integer(kind=kint ), pointer :: INL(:), INU(:)
      integer(kind=kint ), pointer :: IAL(:)
      integer(kind=kint ), pointer :: IAU(:)
  
      real   (kind=kreal):: ALUtmp(3,3)
      integer(kind=kint ):: ip
  
      N = hecMAT%N
      NP = hecMAT%NP
  
negative.code:     integer(kind=kint), dimension(:), allocatable:: inumFI2L, inumFI2U
      integer(kind=kint), dimension(:), allocatable::     FI2L,     FI2U
      real (kind=kreal), dimension(3,3) :: RHS_Aij, DkINV, Aik, Akj
      integer(kind=kint) :: NPLf1,NPLf2,NPUf1,NPUf2,iAS,iconIK,iconKJ
      integer(kind=kint) :: i,jj,ij0,kk,ik,kk1,kk2,L,iSk,iEk,iSj,iEj
      integer(kind=kint) :: icou,icouU,icouU1,icouU2,icouU3,icouL,icouL1,icouL2,icouL3
  