C ZIRO      SOURCE    CHAT      05/01/13    04:23:01     5004
c******************************************************************
c       sous routine utilisee dans l'ecoulement plastique siuvant *
c       la loi de gurson                                          *
c******************************************************************
        subroutine ziro(dphi,wrkgur)
c       this subroutine solve for phi1 the solution of
c       phi=exp(3.|P|./2/sigbar)
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)
*       common /toler/ conv,tol1,tol2
c
c***    init data
        phimax=1.d0-rho/rho0*(1.d0-phi0)
        phi1=phimax
        p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho))
        f=log(phin)+1.5d0*b/sigbar*(1-phin)*(phin-phimax)/(1.d0-phimax)
c
c***    material in compression
        if ( p .gt. 0.d0) then
c
c---    test if there is a root inside the intervall
          if ( f .le. 0.d0) then
            dphi=0.d0
            return
          endif
c
c---    test if the answer will be smaller than phi0
          f=log(phi0)+1.5d0*b/sigbar*(1-phi0)
     &          *(phi0-phimax)/(1.d0-phimax)
          if (f .gt. 0.d0) then
            dphi=phi0*0.99d0-phin
            return
          else
            if (phimax .lt. phi0) phi1=phi0
          endif
c
c---    iterate the tangent approximation
          do 100 iter=1,50
          f=log(phi1)+1.5d0*b/sigbar*(1-phi1)*(phi1-phimax)
     &            /(1.d0-phimax)
c
c---    check for convergence
          syp=sqrt(abs(1.d0+phi1*phi1-2.d0*phi1*cosh(-1.5d0*b*(1-phi1)*
     &         (phi1-phimax)/(1.d0-phimax)/sigbar)))
c
          if ( syp .le. tol1) then
            dphi=phi1-phin
            return
          endif
c---
          fprime=1.d0/phi1+1.5d0*b/sigbar*(1.d0+phimax-2.d0*phi1)/
     &          (1.d0-phimax)
          phi1=phi1-f/fprime
c
  100     continue
c
c***    material in tension
        else
          a1=1+phimax
c
c---    loop over the process
          do 200 iter=1,90
c
          a2=2.d0/3.d0*sigbar/b*(1.d0-phimax)*log(phi1)+phimax
          phi1=(a1-sqrt(a1*a1-4.d0*a2))*0.5d0
c
c---    test if the search is outside the interval
          if ( phi1 .le. phin) then
            dphi=0.d0
            return
          endif
c
c---    check for convergence
          syp=sqrt(abs(1.d0+phi1*phi1-2.d0*phi1*cosh(-1.5d0*b*(1-phi1)*
     &         (phi1-phimax)/(1.d0-phimax)/sigbar)))
c
          if (syp .le. tol1) then
            dphi=phi1-phin
            return
          endif
c
  200     continue
c
        endif
        dphi=phi1-phin
c       the subroutine didnot converge return anyway
        return
c
        end