C SPLIT     SOURCE    PV        08/09/23    21:15:07     6164
c*****************************************************************************
c       sous rouine utilisee dans l'ecoulement plastique de la loi
c       de GURSON
c*****************************************************************************
        subroutine split(dphi,wrkgur)
c
c       this subroutine split up the strain increment in two if a normal
c       evolution reaches a state where phin < phi0
c       xi*de is the first strain increment that makes reaches phin=phi0
c
c       rho and sqrtj2 have been updated in newsta
c
c---    variables
      IMPLICIT INTEGER(I-N)
        IMPLICIT REAL*8 (A-H,O-Z)
        segment wrkgur
         real*8 sigbar, sy0,phi0,rho0,g,b,h
         real*8 epn,phin,sqrtj2,rho,sig(6)
         real*8 e(7),dt
         real*8 conv,tol1,tol2
        endsegment
*        common/prop/sigbar,sy0,phi0,rho0,g,b,h
*        common/state/epn,phin,sqrtj2,rho,sig(6)
c               e is the strain rate:
c                       e(1-6) is the deviatoric part
c                       e(7) is the trace/3
c               dt ist the time increment
*        common /delta/ e(7),dt
*       common /toler/ conv,tol1,tol2
c
c***    set the bounds
        p=(1.d0-phi0)*b*(1.d0-rho0/rho)
        syp=1.d0+phi0*phi0-2.d0*phi0*cosh(1.5d0*p/sigbar)
c
        if ( syp .le. 0.d0) then
c         xi2 is computed such that syp(xi2)=0
          xi2=1.d0+(2.d0/3.d0*sigbar*log(phi0)+(1.d0-phi0)*b*
     &       (1.d0-rho0/rho))
     &       /(1.d0-phi0)/b/e(7)/dt/3.d0
        else
          xi2=1.d0
        endif
        xi1=0.d0
c
        call funk(xi1,f1,wrkgur)
        call funk(xi2,f2,wrkgur)
c
        if ( abs(f2) .le. conv) then
          xi=xi2
          goto 200
        endif
        if ( abs(f1) .le. conv) then
          xi=xi1
          goto 200
        endif
c
        if ( (f2*f1) .gt. 0.d0) then
          write (6,*)  'f2*f1>0 in split'
        endif
c
c***    start iterative process (Ridder's method)
        do 100 j=1,99
c
          xi3=(xi1+xi2)*0.5d0
          call funk (xi3,f3,wrkgur)
c
          dum=xi3+(xi3-xi1)*sign(1.d0,f1-f2)*f3/sqrt(abs(f3**2-f2*f1))
          call funk(dum,fdum,wrkgur)
c
          if (fdum .lt. 0.d0) then
            xi2=dum
            f2=fdum
          else
            xi1=dum
            f1=fdum
          endif
c
          if ( abs(fdum) .le. conv) then
            xi=dum
            goto 200
          endif
c
  100   continue
c
c       split did not converge!'
        xi=xi1
c
  200   continue
c
c***    update new state for xi*de strain increment
        dphi=phi0-phin
        phin=phi0
        p=(1.d0-phi0)*b*(1.d0-rho0/rho)+(1-xi)*
     &     (1.d0-phi0)*b*e(7)*dt*3.d0
        syn=sy0+h*epn
        syp=1.d0+phi0**2-2.d0*phi0*cosh(-1.5d0*p/sigbar)
c
        if (abs(syp) .gt. tol1*tol1) then
          a=1.5d0*phi0*(1.d0-phi0)*sinh(-1.5d0*p/sigbar)
     &          /sqrt(syp)
c
          ratio=syn/(2.d0*h)
          dep=sqrt(ratio**2+abs(dphi/a)*sigbar/h)-ratio
c
          epn=epn+dep
          sy=(syn+h*dep)*sqrt(syp)
          sqrtj2=sy
c
          fac=1.d0+3.d0*g*dep/sqrtj2
c
          do 250 i=1,6
          sig(i)=(sig(i)-(1.d0-xi)*2.d0*g*e(i)*dt)/fac
  250     continue
c
        else
          sqrtj2=0.d0
        endif
c
c***    evolution under (1-xi)*de
        sum=0.d0
c
        do 400 i=1,6
          e(i)=(1.d0-xi)*e(i)
          sig(i)=sig(i)+2*g*e(i)*dt
          sum=sum+sig(i)*sig(i)
          if (i .gt. 3) then
            sum=sum+sig(i)*sig(i)
          endif
  400   continue
        sqrtj2=sqrt(1.5d0*sum)
c
        syn=sy0+h*epn
c
c---    trial state with respect to yield surface
        if ( sqrtj2 .gt. syn )then
          call prandt(wrkgur)
        endif
c
        dphi=1.0d0
        return
        end