anchor.diff:        t(4) = t(4) + zz(1)*Zln(3,k) + zz(4)*Zln(6,k) + zz(7)*Zln(9,k)
         t(5) = t(5) + zz(2)*Zln(3,k) + zz(5)*Zln(6,k) + zz(8)*Zln(9,k)
         t(6) = t(6) + zz(3)*Zln(3,k) + zz(6)*Zln(6,k) + zz(9)*Zln(9,k)
  -    ENDDO
  -    DO l = 1, 6
  +    enddo
  +    do l = 1, 6
         Diag(l,Ic) = Diag(l,Ic) - t(l)
  -    ENDDO
  -    CALL INV3(Diag(1,Ic),ir)
  +    enddo
  +    call INV3(Diag(1,Ic),ir)
       Nch(Ic) = -1
       kk = Par(Ic)
       Nch(kk) = Nch(kk) - 1
  -  END SUBROUTINE S3UM
  +  end subroutine S3UM
   
     !======================================================================!
     !> @brief S3UM1
     !======================================================================!
  -  SUBROUTINE S3UM1(Ic,Xlnzr,Colno,Zln,Temp,Indx)
  -    IMPLICIT NONE
  -    !------
  -    INTEGER, INTENT(IN):: Ic
  -    INTEGER, INTENT(IN):: Xlnzr(*)
  -    INTEGER, INTENT(IN):: Colno(*)
  -    INTEGER, INTENT(OUT):: Indx(*)
  -    REAL(KIND=8), INTENT(OUT):: Temp(9,*)
  -    REAL(KIND=8), INTENT(OUT):: Zln(9,*)
  -    !------
  -    INTEGER:: j
  -    INTEGER:: jc
  -    INTEGER:: jj
  -    INTEGER:: k
  -    INTEGER:: ke
  -    INTEGER:: ks
  -    INTEGER:: l
  -    REAL(KIND=8):: s(9)
  +  subroutine S3UM1(Ic,Xlnzr,Colno,Zln,Temp,Indx)
  +    implicit none
  +    !------
  +    integer, intent(IN):: Ic
  +    integer, intent(IN):: Xlnzr(*)
  +    integer, intent(IN):: Colno(*)
  +    integer, intent(OUT):: Indx(*)
  +    real(KIND=8), intent(OUT):: Temp(9,*)
  +    real(KIND=8), intent(OUT):: Zln(9,*)
  +    !------
  +    integer:: j
  +    integer:: jc
  +    integer:: jj
  +    integer:: k
  +    integer:: ke
  +    integer:: ks
  +    integer:: l
  +    real(KIND=8):: s(9)
   
       ks = Xlnzr(Ic)
       ke = Xlnzr(Ic+1)
       s(1:9) = 0.0D0
  -    DO k = ks, ke - 1
  +    do k = ks, ke - 1
         jc = Colno(k)
         Indx(jc) = Ic
  -      DO jj = Xlnzr(jc), Xlnzr(jc+1) - 1
  +      do jj = Xlnzr(jc), Xlnzr(jc+1) - 1
           j = Colno(jj)
  -        IF ( Indx(j)==Ic ) THEN
  +        if ( Indx(j)==Ic ) then
             s(1) = s(1) + Temp(1,j)*Zln(1,jj) + Temp(4,j)*Zln(4,jj) + Temp(7,j)*Zln(7,jj)
             s(2) = s(2) + Temp(2,j)*Zln(1,jj) + Temp(5,j)*Zln(4,jj) + Temp(8,j)*Zln(7,jj)
             s(3) = s(3) + Temp(3,j)*Zln(1,jj) + Temp(6,j)*Zln(4,jj) + Temp(9,j)*Zln(7,jj)
  
positive.code:         Diag(4,ic) = Diag(4,ic) - Temp(jj,1)*Zln(3,k) - Temp(jj,4)*Zln(6,k) - Temp(jj,7)*Zln(9,k)
          Diag(5,ic) = Diag(5,ic) - Temp(jj,2)*Zln(3,k) - Temp(jj,5)*Zln(6,k) - Temp(jj,8)*Zln(9,k)
          Diag(6,ic) = Diag(6,ic) - Temp(jj,3)*Zln(3,k) - Temp(jj,6)*Zln(6,k) - Temp(jj,9)*Zln(9,k)
        ENDDO
        DO jc = Nstop, ic - 1
          joc = joc + 1
          j1 = Xlnzr(jc)
          j2 = Xlnzr(jc+1)
          DO jj = Xlnzr(jc), Xlnzr(jc+1) - 1
            j = Colno(jj)
            IF ( Indx(j)==ic ) THEN
              Dsln(1,joc) = Dsln(1,joc) - Temp(j,1)*Zln(1,jj) - Temp(j,4)*Zln(4,jj) - Temp(j,7)*Zln(7,jj)
              Dsln(2,joc) = Dsln(2,joc) - Temp(j,2)*Zln(1,jj) - Temp(j,5)*Zln(4,jj) - Temp(j,8)*Zln(7,jj)
              Dsln(3,joc) = Dsln(3,joc) - Temp(j,3)*Zln(1,jj) - Temp(j,6)*Zln(4,jj) - Temp(j,9)*Zln(7,jj)
  
negative.code: 
      !* Solve
      !* Backsubstitute
      call NUSOL0(hecMAT%B,ir)
      call PTIME(t5)
      !* Errors 4
      if ( ir/=0 ) then
        !WINDEBUG
        write (Ifmsg,*) 'error in nusol0. irr = ', ir
        stop
      endif
      call HECMW_MAT_DUMP_SOLUTION(hecMAT)
    end subroutine HECMW_SOLVE_DIRECT
  
    !======================================================================!
    !> @brief PTIME
    !======================================================================!
    subroutine PTIME(Cputim)
      use HECMW_UTIL
      implicit none
      !------
      real(KIND=8), intent(OUT):: Cputim
      !------
      ! cpu time by hour
      Cputim = HECMW_WTIME()
    end subroutine PTIME
  
    !======================================================================!
    !> @brief SETIJ
    !======================================================================!
    subroutine SETIJ(hecMESH,hecMAT)
      use HECMW_UTIL
      implicit none
      !------
      type (HECMWST_LOCAL_MESH), intent(IN)::hecMESH
      type (HECMWST_MATRIX), intent(IN)::hecMAT
      !------
      integer:: i
      integer:: ierr
      integer:: j
      integer:: k
      integer:: kk
      integer:: ntotal
      integer(KIND=kint):: numnp
      integer(KIND=kint):: ndof
      integer(KIND=kint):: ndof2
  
      numnp = hecMAT%NP
      ndof = hecMESH%N_DOF
      ntotal = numnp*ndof
  
      !*NUFACT variables
      NEQns = numnp
      NDEg = ndof
      NTTbr = hecMAT%NP + hecMAT%NPL
      !+hecMAT%NPU if unsymmetric
      ISYm = 0
  
      !*Allocations
      allocate (IROw(NTTbr),STAT=ierr)
      allocate (JCOl(NTTbr),STAT=ierr)
      if ( ierr/=0 ) stop "Allocation error: irow/jcol"
  
      kk = 0
      ndof2 = ndof*ndof
      do j = 1, numnp
        !*Diagonal
        kk = kk + 1
        IROw(kk) = j
        JCOl(kk) = j
        !*Lower
        do k = hecMAT%INDEXL(j-1) + 1, hecMAT%INDEXL(j)
          i = hecMAT%ITEML(k)
          kk = kk + 1
          IROw(kk) = j
          JCOl(kk) = i
        enddo
      enddo
    end subroutine SETIJ
  
    !======================================================================!
    !> @brief MATINI initializes storage for sparse matrix solver.
  