C FUNK      SOURCE    CHAT      05/01/13    00:11:42     5004
c*****************************************************************************
c       sous rouine utilisee dans l'ecoulement plastique de la loi
c       de GURSON
c*********************************************************************
        subroutine funk(xi,f,wrkgur)
c
c       this subroutine compute f for a given xi
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
        dimension s(6)
c
c---    executable

        p=(1.d0-phi0)*b*(1.d0-rho0/rho)+
     &          (1.d0-xi)*(1.d0-phi0)*b*3.d0*e(7)*dt
        syn=sy0+h*epn
        syp=abs(1.d0+phi0**2-2.d0*phi0*cosh(-1.5d0*p/sigbar))
c
c---    plastic strain increment
        if (syp .ge. tol1*tol1 ) then
          a=1.5d0*phi0*(1.d0-phi0)*sinh(-1.5d0*p/sigbar)/sqrt(syp)
          ratio=syn/(2.d0*h)
          dep=sqrt(ratio**2+abs((phi0-phin)*sigbar/a)/h)-ratio
        else
            dep=0.d0
        endif
c
c---    equivalent trial stress
        sum=0.d0
        do 10 i=1,6
                s(i)=sig(i)-2.d0*g*dt*(1.d0-xi)*e(i)
                sum=sum+s(i)**2
                if (i .gt. 3) then
                sum=sum+s(i)**2
                endif
   10   continue
        seqtr=sqrt(1.5d0*sum)
c
        f=(syn+h*dep)*sqrt(syp) + 3.d0*g*dep - seqtr
        f=f/sy0
c
        return
        end
c*******************************************************************************