C SIGF3DSU  SOURCE    CB215821  16/04/21    21:18:25     8920
      SUBROUTINE sigf3d_sub(esigmaesave,x1,x2,sigmaf,H66,P2,pi2,
     .              ag,i6,fc,fb,parahot3,idimpara3,lerror,lcp,
     .              esigmae1,U,ldivide,n,xsub1,xsub2,esigmae)

c     This subroutine deals with the subdivision of the plastic increments for
c       improving the convergence of the subroutine sigf3d_implicit (algorithm
c       to find the stress sigmaf from the trial stress esigmae and the plastic
c       increments x1 and x2)

      IMPLICIT REAL*8 (A-B,D-H,O-Z)
      implicit integer (I-K,M,N)
      implicit logical (L)
      implicit character*10 (C)

      dimension sigmaf(6),sigmai(6)
c       final stress vector
      dimension esigmae(6),esigmae1(6),dsigma(6)
c       trial stress vector
      dimension H66(6,6),pi2(6)
      dimension dgc(6),dgt(6),vloc(6),trav6(6),P2(6,6)
      dimension parahot3(idimpara3),rcossigmapr(3,3),B(6,6)
      dimension d2gc(6,6),d2gt(6,6),trav66(6,6),h(6),U(6,6)
      dimension esigmaesave(6)

      i0 = 0
      i1 = 1
      i2 = 2
      i3 = 3
      i4 = 4
      i5 = 5
      i7 = 7

      r2 = 2.

      ldivide = .false.
      Dlambdat = xsub1
      Dlambdac = xsub2
      do jloc=1,6
        do iloc=1,6
          U(iloc,jloc) = 0.d0
          d2gc(iloc,jloc) = 0.d0
          d2gt(iloc,jloc) = 0.d0
        enddo
      enddo
      do iloc=1,6
        U(iloc,iloc) = 1.0d0
      enddo

c     Case Dlambdat = 0 and Dlambdac = 0
      if ((Dlambdat.eq.0.).and.(Dlambdac.eq.0.)) then
        do iloc=1,6
          sigmaf(iloc) = esigmae(iloc)
        enddo
        return
      endif

c     INITIALISATION
c     increment of stress
      do iloc=1,6
        dsigma(iloc) = 0.d0
      enddo
c     effective stress for the first iteration
      do iloc=i1,i6
        sigmai(iloc) = esigmae(iloc)
      end do
      if (lcp) then
        sigmai(3) = 0.d0
        sigmai(5) = 0.d0
        sigmai(6) = 0.d0
      endif
c     function h
c     case first iteration (sigmai = 0)
      rloc = 0.d0
      do iloc=1,6
        rloc = rloc + ABS(sigmai(iloc))
      enddo
c     if (rloc.lt.10.) then
      if (rloc.eq.0.) then
        call dsiggc(esigmae,dgc,ag,P2,pi2,i1,i6,fc,fb,lcp)
        call dsigRank(esigmae,dgt,i3,i6,parahot3,idimpara3,
     .                rkappait,lcp)
c       effective stress
        do iloc=1,6
          sigmaf(iloc) = esigmae(iloc)
        end do
      else
        call dsiggc(sigmai,dgc,ag,P2,pi2,i1,i6,fc,fb,lcp)
        call dsigRank(sigmai,dgt,i3,i6,parahot3,idimpara3,
     .                rkappait,lcp)
c       effective stress
        do iloc=1,6
          sigmaf(iloc) = sigmai(iloc)
        end do
      endif
      do iloc=i1,i6
        vloc(iloc) = Dlambdat*dgt(iloc) + Dlambdac*dgc(iloc)
      end do
      call mulAB(H66,VLOC,trav6,i6,i6,i6,i1)
      do iloc=1,6
        h(iloc) = sigmaf(iloc) - esigmae(iloc) + trav6(iloc)
      enddo
c     Jacobien
      if (Dlambdac.ne.0)
     .   call ddsiggc(sigmaf,d2gc,ag,P2,pi2,i1,i6,fc,fb,lcp)

      if (Dlambdat.ne.0)
     .   call ddsigrank(sigmaf,d2gt,i1,i2,i3,i4,i5,i6,lcp,
     .                  parahot3,idimpara3)
      do jloc=1,6
        do iloc=1,6
          B(iloc,jloc) = Dlambdat * d2gt(iloc,jloc) +
     .                   Dlambdac * d2gc(iloc,jloc)
        enddo
      enddo
      call mulAB(H66,B,trav66,i6,i6,i6,i6)
      do jloc=1,6
        do iloc=1,6
          B(iloc,jloc) = U(iloc,jloc) + trav66(iloc,jloc)
        enddo
      enddo

c     counter
      iter2 = 1

c     ITERATIVE SOLVER
c     function to solve
  15  continue
      rloc = 0.
      do iloc=1,6
        rloc = rloc + ABS(h(iloc))
      enddo
      if ((iter2.eq.1).or.((rloc.gt.10.).and.(iter2.lt.30))) then

        call invert9(B,6,6)
        do jloc=1,6
          do iloc=1,6
            trav66(iloc,jloc) = -B(iloc,jloc)
          enddo
        enddo

        call mulAB(trav66,h,dsigma,6,6,6,1)
        do iloc=1,6
          sigmaf(iloc) = sigmaf(iloc) + dsigma(iloc)
        enddo
c       calculation of dgt and dgc
c       test on sigma - if sigma close to 0, use sigma,tr for calculation of dgt and dgc
        rloc = 0.d0
        if (lcp) then
          rloc = (sigmaf(1))**2 + (sigmaf(2))**2 + (sigmaf(4))**2
        else
          rloc = (sigmaf(1))**2 + (sigmaf(2))**2 + (sigmaf(3))**2 +
     .           (sigmaf(4))**2 + (sigmaf(5))**2 + (sigmaf(6))**2
        endif
        if (rloc.gt.0.) then
          call dsiggc(sigmaf,dgc,ag,P2,pi2,i1,i6,fc,fb,lcp)
          call dsigRank(sigmaf,dgt,i3,i6,parahot3,idimpara3,
     .                  rkappait,lcp)
        else
          call dsiggc(esigmae,dgc,ag,P2,pi2,i1,i6,fc,fb,lcp)
          call dsigRank(esigmae,dgt,i3,i6,parahot3,idimpara3,
     .                  rkappait,lcp)
        endif
c       function h
        do iloc=i1,i6
          vloc(iloc) = Dlambdat*dgt(iloc) + Dlambdac*dgc(iloc)
        end do
        call mulAB(H66,VLOC,trav6,i6,i6,i6,i1)
        do iloc=1,6
          h(iloc) = sigmaf(iloc) - esigmae(iloc) + trav6(iloc)
        enddo
c       Jacobien
        if (Dlambdac.ne.0)
     .     call ddsiggc(sigmaf,d2gc,ag,P2,pi2,i1,i6,fc,fb,lcp)
        if (Dlambdat.ne.0)
     .     call ddsigrank(sigmaf,d2gt,i1,i2,i3,i4,i5,i6,lcp,
     .                    parahot3,idimpara3)
        do jloc=1,6
          do iloc=1,6
            B(iloc,jloc) = Dlambdat*d2gt(iloc,jloc) +
     .                     Dlambdac * d2gc(iloc,jloc)
          enddo
        enddo
        call mulAB(H66,B,trav66,i6,i6,i6,i6)
        do jloc=1,6
          do iloc=1,6
            B(iloc,jloc) = U(iloc,jloc) + trav66(iloc,jloc)
          enddo
        enddo
        iter2 = iter2+1

        go to 15
      endif

      if (iter2.eq.30) then
        ldivide = .true.
      endif

      RETURN
      END