C NEWSTA    SOURCE    PV        15/04/13    21:15:15     8474
c*****************************************************************************
c       sous rouine utilisee dans l'ecoulement plastique de la loi
c       de GURSON
c**********************************************************************
        subroutine newsta(wrkgur)
c       this subroutine update the material state given the strain increment e
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---    executable
c
c       parametres pour la convergence et les tolerances d'erreur
        conv=1.d-4
        tol1=1.d-5
        tol2=1.d-6
c
c---    trial stress
        dum=(rho0-rho)/rho+e(7)*3.d0*dt+1.d0
        rho=rho0/dum
        sum=0.d0
        do 10 i=1,6
             sig(i)=sig(i)+2.d0*g*e(i)*dt
             sum=sum+sig(i)*sig(i)
             if (i .gt. 3) then
               sum=sum+sig(i)*sig(i)
             endif
   10   continue
        sqrtj2=sqrt(1.5d0*sum)
        p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho))
        syn=sy0+h*epn
c
c***    test if phin < phi0 a Prandtl-Reuss model is used
        if ((phin .le. phi0) .and. ( p .ge. 0.d0)) then
          if (sqrtj2 .gt. syn) then
            call prandt(wrkgur)
          endif
          return
        endif
c
        syach=-1.5d0*p/sigbar
*  borner pour eviter un depassement
        syach=min(675d0,max(-675d0,syach))
        sy=syn**2*(1.d0+phin**2-2.d0*phin*cosh(syach))
        yield=sqrtj2*sqrtj2-sy
C
c---    initialize paramter
        dphi=0.d0
        iter=0
c
c***    call solve if plastic evolution
        if (yield .gt. 0.d0) then
          call solve(dphi,iter,wrkgur)
c
c***    check if increment has been splitted in zero()
c               this implies phin=phi0
          if ( dphi .ge. 1.d0)  then
            p=(1.d0-phi0)*b*(1.d0-rho0/rho)
            goto 60
          endif
c
c---    update new state value given new dphi
          syn=sy0+h*epn
          phin=phin+dphi
          p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho))
          syp=1.d0+phin**2-2.d0*phin*cosh(-1.5d0*p/sigbar)
          dep=(sqrtj2-syn*sqrt(syp))/(3*g+h*sqrt(syp))
c
          if (abs(syp) .ge. tol1*tol1) then
c           a=1.5d0*phin*(1.d0-phin)*sinh(-1.5d0*p/sigbar)/sqrt(syp)
c
c           if ( a .le. tol2) then
c               dep=(sqrtj2-syn*(1-phin))/(3*g+h*(1-phin))
c           else
c                ratio=syn/(2.d0*h)
c                dep=sqrt(ratio*ratio+abs(dphi/a)*sigbar/h)-ratio
c           endif
c
            epn=epn+dep
            sy=(syn+h*dep)*sqrt(syp)
            sqrtj2=sy
c
c---        return stress to yield surface
            fac=1.d0+3.d0*g*dep/sqrtj2
            do 20 i=1,6
              sig(i)=sig(i)/fac
   20       continue
          else
            sqrtj2=0.d0
          endif
        endif
c
  60    continue
c
c***    echo new state
        return
c
        end