C FUNC      SOURCE    CHAT      05/01/13    00:11:36     5004
c*****************************************************************************
c       sous rouine utilisee dans l'ecoulement plastique de la loi
c       de GURSON
c**********************************************************************
        subroutine func(f,q,wrkgur)
c       this subroutine compute the value of the function
c       which zero is the increment of phi
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
*       common /toler/ conv,tol1,tol2
c
c***    evaluate function -- q is dphi
        phi=phin+q
        syn=sy0+h*epn
        p=(1.d0-phi)*b*(1.d0-(1.d0-phi)*rho0/((1.d0-phi0)*rho))
        syp=abs(1.d0+phi**2-2.d0*phi*cosh(-1.5d0*p/sigbar))
c
        if (syp .ge. tol1*tol1 ) then
        a=1.5d0*phi*(1.d0-phi)*sinh(-1.5d0*p/sigbar)/sqrt(syp)
c
            if ( abs(a) .le. (1.d-2 * tol2) ) then
c       this is not consistent but this speeds up the convergence
                dep=(sqrtj2-syn*(1-phi))/(3*g+h*(1-phi))
            else
                ratio=syn/(2.d0*h)
                dep=sqrt(ratio**2+abs(q/a)*sigbar/h)-ratio
            endif
c
        else
            dep=0.d0
        endif
c
        syn1p=(syn+h*dep)*sqrt(syp)
        f=syn1p+3.d0*g*dep-sqrtj2
        f=f/sy0
c
        return
        end
c*****************************************************************************