C SOLVE     SOURCE    CB215821  26/05/04    21:15:05     12531          
c*****************************************************************************
c       sous rouine utilisee dans l'ecoulement plastique de la loi
c       de GURSON
c**********************************************************************
        subroutine solve(dphik,iter,wrkgur)
c       this subroutine solve the non linear equation which solution is
c       the phi increment
c
c---    variables
        IMPLICIT INTEGER(I-N)
        IMPLICIT REAL*8 (A-H,O-Z)
        
        
C       Une SUBROUTINE nommee "split" est entree dans la norme FORTRAN 2023
C       GCC ne veut plus compiler cette source car le nombre d'argument de SPLIT
C       ne correspond pas au modele INTRINSIC ==> Ajout de EXTERNAL SPLIT
        EXTERNAL SPLIT
        
        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***    initial data
        p=(1.d0-phin)*b*(1.d0-(1.d0-phin)*rho0/((1.d0-phi0)*rho))
        phimax=1.d0-rho/rho0*(1.d0-phi0)
        if ( phimax .lt. 0.d0) phimax=0.d0
c
c***    find the first dphi for which syp>0
        call ziro(dphik,wrkgur)
        phik=phin+dphik
c
        if ( phik .lt. phi0) then
          call split(dphik,wrkgur)
          return
        endif
c
c---    initialize the search
        call func(fk,dphik,wrkgur)
        iter=1
c
        if ( abs(fk) .le. conv) return
c
        if ( fk .gt. 0.d0) then
           return
        endif
c
c*******perform iterations until find fk>0
        if ( phimax .gt. 0.d0) then
c
c---    compute the constant used to determine dphi
        c=sqrt(abs(phimax-phin)/phimax/h/b)*2.d0*sigbar/3.d0*
     &        (h*(1-phimax)+3.*g)
        c=c/sy0/2.d0
c---    iterates until converges
        do 100 i=1,50
c
c       check for convergence
          if (abs(fk).le. conv) return
c
c       check if fk>0 and then exit loop
          if (fk .gt. 0.d0) goto 200
c
c       check if split is needed
          if ( phik .lt. phi0) then
            call split(dphik,wrkgur)
            return
          endif
c
c---    keep track of the previous step
                fj=fk
                dphij=dphik
c
                iter=iter+1
c
c---    find new value of dphik using the interpolation function
                d=fk-c*sqrt(1.d0/abs(phimax-phik))
                phik=phimax+sign(1.d0,p)*(c/d)**2
                dphik=phik-phin
                call func(fk,dphik,wrkgur)
c
  100   continue
c
c*******case when phimax is negative
        else
c
c---    initialize data for interpolation function
        p0=b*(1-rho0/rho/(1-phi0))
        c2=(h+3.d0*g)*sqrt(phin*sigbar/1.5d0/h
     &      /sinh(1.5d0*p0/sigbar))
     &      /1.d1
c
c---    iterates til converges
        do 150 i=1,50
c
c       check for convergence
          if (abs(fk).le. conv) return
c
c       check if fk>0 and then exit loop
          if (fk .gt. 0.d0) goto 200
c
c       check if split is needed
          if ( phik .lt. phi0) then
            call split(dphik,wrkgur)
            return
          endif
c
c---    keep track of the previous step
          fj=fk
          dphij=dphik
          iter=iter+1
c
c---    find new value of dphik using the interpolation function
          d=fk-c2*sqrt(1.d0/phik)
          phik=sign(1.d0,p)*(c2/d)**2
          dphik=phik-phin
          call func(fk,dphik,wrkgur)
c
  150   continue
        endif
c
c---    the previous algorithm didn't converge
        return
  200   continue
c
c*******find zero using Ridders method
c       see Numerical recipee in Fortran p351
        do 300 i=1,200
           iter=iter+1
c
           half=(dphik+dphij)/2
           call func(fm,half,wrkgur)
           s=sqrt(abs(fm**2-fk*fj))
           dum=half+(half-dphij)*(sign(1.d0,fj-fk)*fm/s)
c
c---    new interval
           call func(fdum,dum,wrkgur)
           if (fdum .lt. 0.d0) then
              dphij=dum
              fj=fdum
           else
              dphik=dum
              fk=fdum
           endif
c
           if ( phin+dphij .lt. phi0) then
             call split(dphik,wrkgur)
             return
           endif
c
c---    check for convergence
        if (abs(fdum) .lt. conv) then
              dphik=dum
              if ( (phin+dphik) .le. phi0) then
                 call split(dphik,wrkgur)
                 return
              endif
              return
        endif
c
  300   continue
        return
c
        end
c*************************************************************