C COUL2     SOURCE    AM        08/08/29    21:15:05     6145
      SUBROUTINE COUL2(IB,IGAU,NSTRS,SIG0,EPIN0,VAR0,NVARI,
     .     DEPST,IFOUR,XMAT,NMATT,IVAL,DD,SIGF,DEFP,VARF,KERRE)
C--------------------------------------------------------------
c
C
C         PLASTICITE MODELE 2D MOHR COULOMB ELEMENTS JOINTS
C
C ENTREES
C     IB              = NUMERO DE L'ELEMENT
C     IGAUSS          = NUMERO DU POINT DE GAUSS
C     NSTRS           = NOMBRE DE COMPOSANTES DE CONTRAINTES
C     SIG0            = CONTRAINTES INITIALES (AU PAS PRECEDENT)
C     EPIN0           = DEFORMATIONS INITIALES (AU PAS PRECEDENT)
C     VAR0(NVARI)     = VARIABLES INTERNES DEBUT
C      VAR0(1)        = DEFORMATION PLASTIQUE EQUIVALENTE
C      VAR0(2)        = EPOUN (DEFO. PLASTIQUE DUE A L'OUVERTURE SEULE)
C      VAR0(3)        = STAT (ETAT DU JOINT)
C      VAR0(4)        = LAM1  ( MULTIPLICATEUR PLASTIQUE EN TRACTION )
C     DSIGT(NSTRS)    = INCREMENT DE CONTRAINTES TOTALES
C     XMAT(NCOMAT)    = COMPOSANTES DE MATERIAU
C     NMATT           = NOMBRE DE COMPOSANTES DE MATERIAU
C     IVAL(NCOMAT)    = INDICE DES COMPOSANTES DE MATERIAU
C     DD(NSTRS,NSTRS) = MATRICE DE HOOK
C
C  SORTIES
C     SIGF(NSTRS)     = CONTRAINTES FINALES
C     DEFP(NSTRS)     = INCREMENTS DE DEFORMATIONS INELASTIQUES FINALES
C     VARF(NVARI)     = VARIABLES INTERNES FINALES
C     KERRE           = INDICE D'ERREUR
C                     = 60 PAS DE CONVERGENCE APRES 10 ITERATIONS
C                          INTERNES (VARIABLE ISSU DANS COULM)
C                     = 61
C
C-----------------------------------------------------------------------
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C    Include contenant quelques constantes dont XPI :
-INC CCREEL
C
      DIMENSION SIG0(*),EPIN0(*),VAR0(*),DEPST(*),XMAT(*)
      DIMENSION DD(NSTRS,NSTRS)
      DIMENSION SIGF(*),DEFP(*),VARF(*)
      KERRE = 0
C
C
C
C ----------PARAMETRES POUR LE TEST DE NULLITE DES CONTRAINTES ET
C ----------DEFORMATIONS
C
      PRESP = +1.D-5
      PRESM = -1.D-5
ccc       PREEP = +1.D-10
ccc       PREEM = -1.D-10
      PREEP = +1.D-4
      PREEM = -1.D-7
c
C
C ----------PARAMETRE D'IMPPRESSION
C
      IMPRI = 0
C
C ----------NOMBRE MAXIMUM D'ITERATIONS INTERNES AUTORISEES
C
      MAXITE = 1000
C
C
C------CHARGEMENT DU EPSILON DE PRECISION DES ITERATIONS INTERNES
C
      PRECIS =  1.D-6
      PRECIN = -1.D-6
C
C
C------CHARGEMENT DE LA MATRICE DE HOOK INVERSE
C------POUR LE CALCUL DE L'INCREMENT DE DEFORMATION
C------A PARTIR DE L'INCREMENT DE CONTRAINTE
C
      DD(1,1) = 1.0D0/XMAT(1)
      DD(1,2) = 0.0D0
      DD(2,1) = 0.0D0
      DD(2,2) = 1.0D0/XMAT(2)
c
c           contrainte en fin de pente ec
c
      sigsei = -xmat(6)*xmat(2)
c
C
C-----------INCREMENT DE DEFORMATION
C
      DEPS1 = DEPST(1)
      DEPSN = DEPST(2)
C
C             TRAITEMENT DU PASSAGE DIRECT ENTRE UN ETAT OUVERT
C             (I.E. EXISTENCE D'UN JEU AU TEMPS N) A UN ETAT
C             COMPRIME AU TEMPS N+1
C
      DEFPN=0.0D0
      INDIC=0
      IFIN =0
      STOTAN = 0.0D0
      DEFP11 = 0.0D0
      DEFP12 = 0.0D0
      ALP    = 1.0D0
      S1INT  = 0.0D0
      ZJEU = EPIN0(2)
C                                CAS DU RAPPROCHEMENT DES LEVRES
C
C  AM  CORRECTION POUR JOINT AVEC RESISTANCE EN TRACTION
C
C        IF (DEPSN.LT.-1.D-4) THEN
      IF (DEPSN.LT.-1.D-4.AND.SIG0(2).GT.PRESP) THEN
C                                IL EXISTE UN JEU
         IF (ZJEU.GT.1.D-3) THEN
C                                PASSAGE DIRECT D'UN ETAT OUVERT
C                                A UN ETAT COMPRIME
            valdif=(ABS(DEPSN) - ZJEU)/ABS(DEPSN)
            IF (valdif.gt.1.D-3) THEN
C
               ALP = ZJEU/(ABS(DEPSN))
               STOTAN= DEPS1
               DEPS1 = ALP * DEPS1
               STOCK = DEPSN
               DEPSN = ALP * DEPSN
               INDIC = 1
               GO TO 99
            ENDIF
C
         ENDIF
C
      ENDIF
C
 98   CONTINUE
C
C                                IL RESTE A PARCOURIR LE TRAJET
C                                (1-ALP) * DEPSN
C
      IF (INDIC.EQ.1) THEN
         DEPS1=(1.0D0-ALP)*STOTAN
         DEPSN=(1.0D0-ALP)*STOCK
         INDIC=0
         IFIN =1
      ENDIF
C
 99   CONTINUE
C
C
C ----CAS GENERAL SANS TRANSITION DIRECTE D'UN ETAT OUVERT A UN ETAT
C ----A UN ETAT COMPRIME
C
C
C
C
      IF ((INDIC.EQ.1).OR.(IFIN.EQ.0)) THEN
C
C--------------------CALCUL DE L ETAT DE CONTRAINTES, DE DEFORMATIONS
C--------------------DES VARIABLES INTERNES AU PAS PRECEDENT
C
C------CHARGEMENT DES VALEURS AU PAS PRECEDENT
C
C
C----------DEFORMATIONS
C
C             DEFORMATIONS INELASTIQUES
C
         EPSS1 = EPIN0(1)
         EPSN1 = EPIN0(2)
C
C             DEFORMATIONS ELASTIQUES
C
         DEFS10 = DD(1,1)*SIG0(1)
c
         if   (sig0(2).lt.sigsei) then
            DEFN0 = ( -xmat(6) ) + ( (sig0(2) - sigsei) / xmat(5) )
         else if ( (sig0(2).ge.sigsei) .and. (sig0(2).le.presp) ) then
            DEFN0 = DD(2,2)*SIG0(2)
         else if   (sig0(2).gt.presp)  then
C
C  AM  CORRECTION POUR JOINT AVEC RESISTANCE EN TRACTION
C           DEFN0 = 0.0D0
            DEFN0 = DD(2,2)*SIG0(2)
         endif
C
C             DEFORMATIONS TOTALES AU PAS PRECEDENT N
C
         EPS1N  = EPSS1 + DEFS10
         EPSNN  = EPSN1 + DEFN0
C
      ENDIF
C
C
C
C
C ----CAS DU 2ND PASSAGE : D'UN ETAT FERME ON PASSE A UN ETAT COMPRIME
C
C
C
C
      IF (IFIN.EQ.1) THEN
C
C--------------------CALCUL DE L ETAT DE CONTRAINTES, DE DEFORMATIONS
C--------------------DES VARIABLES INTERNES AU PAS PRECEDENT
C
C------CHARGEMENT DES VALEURS AU PAS PRECEDENT
C
C
C----------DEFORMATIONS
C
C             DEFORMATIONS INELASTIQUES
C
         EPSS1 = DEFP11
         EPSN1 = 0.0D0
C
C             DEFORMATIONS ELASTIQUES
C
         DEFS10 = DD(1,1)*S1INT
         DEFN0  = 0.0D0
C
C             DEFORMATIONS TOTALES AU PAS PRECEDENT N
C
         EPS1N  = EPSS1 + DEFS10
         EPSNN  = EPSN1 + DEFN0
C
      ENDIF
C
C             ESTIMATION DE LA DEFO. TOTALE AU PAS N+1
C
      EP1NP1 = EPS1N + DEPS1
      EPNNP1 = EPSNN + DEPSN
C
C----------EPSILON ELASTIQUE (1 2 N) AU PAS N+1 A L'ITERATION 0
C
      EE1P10 = DEFS10 + DEPS1
      EENP10 = DEFN0  + DEPSN
C
C----------SIGMA (1 2 N) AU PAS N+1 A L'ITERATION 0
C
      S1NP10 = XMAT(1) * EE1P10
c
      x6 = -1.0d0 * xmat(6)
c
c  AM  CORRECTION POUR JOINT AVEC RESISTANCE EN TRACTION
c
c        if  (depsn.lt.preem) then
      if  (depsn.lt.preem.and.sig0(2).le.presp) then
c
c -----------fermeture du joint
c
         if (epnnp1.le.preep) then
c
c --------------fermeture du joint avec compression des levres
c
            if ( (eenp10.ge.x6) .and. (eenp10.le.preep) ) then
c
c ----------------compression : pente EC
c
               SNNP10 = XMAT(2) * EENP10
            endif
c
            if   (eenp10.lt.x6) then
c
c ----------------compression : pente Ef
c
               snnp10= ( xmat(2) * x6 ) +
     .              ( (eenp10 + xmat(6)) * xmat(5) )
            endif
c
         endif
c
         if ((epnnp1.gt.preem).or.(indic.eq.1)) then
c
c --------------rapprochement des levres sans compression
c --------------ou fermeture des levres sans compression
c
            SNNP10 = 0.0D0
c
         endif
c

c
      else
c
c -----------ouverture du joint ou fermeture si joint pas casse
c
         SNNP10 = XMAT(2) * EENP10
c
      endif
C
C-----------VARIABLES INTERNES
C
C          DEFORMATION PLASTIQUE EQUIVALENTE
C
      EPSEQ=VAR0(1)
C
C          DEFORMATION PLASTIQUE DUE A L'OUVERTURE SEULE
C
      EPOUN=VAR0(2)
C
C          ETAT DU JOINT
C
C            STAT   = 0 LEVRES EN CONTACT, AVEC OU SANS CISAILLEMENT
C            STAT   = 1 JOINT EN OUVERTURE PENTE EC
C            STAT   = 2 JOINT EN FERMETURE PENTE EC
C            STAT   = 3 JOINT EN FERMETURE PENTE EF
C
C
C
      STAT =VAR0(3)
C
C           MULTIPLICATEUR PLASTIQUE EN TRACTION
C
      XLAM1=VAR0(4)
C
C
C----    MOHR-COULOMB NORMAL
C
      C     = XMAT(7)
      PHI   = XMAT(8)
      FTRC  = XMAT(4)
C
C  AM 28/7/95
C  FTRC EST EN 4 ET NON EN 9 A CAUSE DE RHO ET ALPH
C  VOIR COMMENTAIRE DANS ECOUL2
C
C   LA LIMITE EN TRACTION EST MISE BRUTALEMENT A 0
C   DES QU'ELLE EST DEPASSEE.
C
C
      PHI   = (PHI*XPI)/180.D0
ccc       ZMU   = XMAT(8)
      ZMU   = 0.0D0
      ZMU   = (ZMU*XPI)/180.D0
C
C     MISE A 0 DE LA LIMITE EN TRACTION SI DEJA CASSE OU SI JEU
C
      RTRAC=FTRC
      IF(XLAM1.GT.0.D0) RTRAC=0.D0
C
C-----------ETUDE DE L'ECOULEMENT SELON LES DIFFERENTS CRITERES
C
      S2 = 0.0D0
      CRICIS = CALG(S1NP10,S2,SNNP10,PHI,C)
      CRINOR = SNNP10 - RTRAC
      CRICIP = CALN(S1NP10,S2,SNNP10,PHI,C)
      ABSTAU = SQRT( (S1NP10 ** 2) + (S2 ** 2) )
C
C----------AIGUILLAGE SUIVANT LA ZONE 0, 1, 2, 3
C
      IF ((CRICIS.LE.PRECIS).AND.(CRINOR.LE.PRECIS)) THEN
         GO TO 10
      ENDIF
C
      IF ((CRICIS.GT.PRECIS).AND.(CRICIP.LE.PRECIS)) THEN
         GO TO 40
      ENDIF
C
      IF ((CRICIP.GT.PRECIS).AND.(ABSTAU.GT.C)) THEN
         GO TO 20
      ENDIF
C
      IF ((ABSTAU.LE.C).AND.(CRINOR.GT.PRECIN)) THEN
         GO TO 30
      ENDIF
C
C
C-----------CAS DE LA CONVERGENCE A L ITERATION 0 :
C
 10   CONTINUE
C
C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
C
C
      IF (IMPRI.EQ.1) THEN
         WRITE(6,*) 'PASSAGE N. 1'
         WRITE(6,*) 'LES DEUX CRITERES SONT VERIFIES-PAS D ECOULEMENT'
         WRITE(6,*) 'CONVERGENCE A L ITERATION 0'
      ENDIF
C
C
C
      IF (INDIC.EQ.1) THEN
C                              FERMETURE DU JOINT SANS COMPRESSION
         S1INT = S1NP10
         DEFP11= 0.0D0
cc             DEFPN = 0.0D0
         DEFPN = -epin0(2)
         GO TO 98
C
      ELSE
C
         IF (IFIN.EQ.0) THEN
C
            DEFP(1) = 0.0D0
C
            if ( abs(snnp10).le.presp) then
c
c            deformation plastique du pas (n+1) moins
c            deformation plastique du pas (n)
c
               defp(2) = epnnp1 - epsn1
            else
c
c            le joint se deforme toujours elastiquement en fermeture
c
               DEFP(2) = 0.0D0
            endif
c
         ELSE
c
c             cas du passage direct d'un etat ouvert a un etat
c             comprime
c
            DEFP12  = 0.0D0
            DEFP(1) = 0.0d0
cc                DEFP(2) = 0.0d0
            DEFP(2) = defpn
            IFIN = 0
         ENDIF
C
C     CALCUL DES VARIABLES INTERNES : DEPEQ  STAT  EPOUN
C
c
c      calcul du delta epsilon plastique entre
c      le pas actuel et le pas precedent
c
         delep1 = defp(1)
         delepn = defp(2)
c
c      produit contracte pour le calcul de la defo. plasti. cumulee
c
         pro1 = ( delep1 * delep1 ) + ( delepn * delepn )
         depeq = sqrt( (2.0D0/3.0D0) * pro1 )
         VARF(1) = EPSEQ + DEPEQ
         VARF(4) = XLAM1
c
c      calcul de la deformation normale due a l'ouverture seule
c
         if ( abs(snnp10).le.presp) then
            epoun = defp(2)
            VARF(2) = EPOUN
         else
            epoun = 0.0d0
            VARF(2) = EPOUN
         endif
c
c     etat du joint     stat = 1   cisaillement pur, levres en contact
c                               ou ouverture du joint pente EC
c                       stat = 2   compression du joint pente EC
c                       stat = 3   compression du joint pente EF
c
         if ( abs(snnp10).le.presp) then
            stat = 1.D0
            if ( abs(epnnp1).le.1.D-8 ) then
               stat = 0.D0
            endif
            VARF(3) = STAT
         endif
         if ( (snnp10.gt.sigsei).and.(snnp10.lt.presm) ) then
            stat = 2.D0
            VARF(3) = STAT
         endif
         if (snnp10.le.sigsei) then
            stat = 3.D0
            VARF(3) = STAT
         endif
C
         SIGF(1) = S1NP10
         SIGF(2) = SNNP10
C
         RETURN
c
      ENDIF
C
C
C----------CRITERES DE L'EFFORT NORMAL ET DE CISAILLEMENT
C----------TOUS LES DEUX VIOLES EN MEME TEMPS
C
 20   CONTINUE
C
C            ECOULEMENT SELON LES DEUX CRITERES A L'ITERATION ITER
C
      ITER = 1
      DELAM1=0.D0
C
 2    CONTINUE
C
      IF (ITER.GT.MAXITE) THEN
         WRITE(6,*) ' LES DEUX CRITERES SONT VIOLES '
         WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES'
         KERRE = 2
         RETURN
      ENDIF
C
      IF (IMPRI.EQ.1) THEN
         write(6,*) 'PASSAGE N. 2'
         write(6,*) 'LES DEUX CRITERES SONT VIOLES'
         write(6,*) 'ITERATION',ITER
      ENDIF
C
C
C            CALCUL DES CORRECTEURS PLASTIQUES
C
C
C                DERIVEE DU CRITERE F SUIVANT SIGMA
C
      DFIDS1 = 0.0D0
      DFIDSN = 1.0D0
C
C                DERIVEE DU CRITERE G SUIVANT SIGMA
C
      if ( s1np10.lt.precis ) then
         dgids1 = -1.d0
      else
         dgids1 =  1.d0
      endif
c
      DZGIDS = TAN(ZMU)
      DGIDSN = TAN(PHI)
C
C                MATRICE DE HOOK * DFI/DSIGMA
C
      HODF1 = XMAT(1) * DFIDS1
      HODFN = XMAT(2) * DFIDSN
C
C                MATRICE DE HOOK * DGI/DSIGMA
C
      HODG1 = XMAT(1) * DGIDS1
      HODGN = XMAT(2) * DZGIDS
C
C                FIFI = DFI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
C
      FIFI = (DFIDS1 * HODF1) + (DFIDSN * HODFN)
C
C                FIGI = DFI/DSIGMA * MATRICE DE HOOK * DGI/DSIGMA
C
      FIGI = (DFIDS1 * HODG1) + (DFIDSN * HODGN)
C
C                GIFI = DGI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
C
      GIFI = (DGIDS1 * HODF1) + (DGIDSN * HODFN)
C
C                GIFI = DGI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
C
      GIGI = (DGIDS1 * HODG1) + (DGIDSN * HODGN)
C
C                RESOLUTION DU SYSTEME POUR TROUVER LES CORRECTEURS
C                PLASTIQUES DLAM1 ET DLAM2
C
      DET = (FIFI * GIGI) - (FIGI * GIFI)
C
      IF (ABS(DET).LE.PRECIS) THEN
         WRITE(6,*) ' SYSTEME INDETERMINE (F G): ITERATION ',ITER
         KERRE = 2
         RETURN
      ENDIF
C
      DLAM1 = (CRINOR * GIGI) - (CRICIS * FIGI)
      DLAM1 = DLAM1 / DET
      DELAM1 = DELAM1 + DLAM1
      RTRAC=0.D0
C
      DLAM2 = (CRICIS * FIFI) - (CRINOR * GIFI)
      DLAM2 = DLAM2 / DET
C
C                CALCUL DU NOUVEL ETAT DE CONTRAINTES
C
      S1NP11=S1NP10 - (XMAT(1) * ( (DLAM1*DFIDS1) + (DLAM2*DGIDS1)))
      SNNP11=SNNP10 - (XMAT(2) * ( (DLAM1*DFIDSN) + (DLAM2*DZGIDS)))
C
C    TEST DE LA CONVERGENCE
C
      S2 = 0.0D0
      CRICIS = CALG(S1NP11,S2,SNNP11,PHI,C)
      CRINOR = SNNP11 - RTRAC
C
      IF ((ABS(CRICIS).LE.PRECIS).AND.(ABS(CRINOR).LE.PRECIS)) THEN
C
C----------LES DEUX CRITERES SONT VERIFIES
C
C
         IF (IMPRI.EQ.1) THEN
            write(6,*) 'LES DEUX CRITERES SONT VERIFIES'
            write(6,*) 'CONVERGENCE A L ITERATION',ITER
         ENDIF
C
C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
C
C
         SIGF(1) = S1NP11
         SIGF(2) = SNNP11
C
         DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1
c
         if ( abs(sigf(2)) .le. presp ) then
c
c            deformation plastique du pas (n+1) moins
c            deformation plastique du pas (n)
c
            DEFP(2) = EPNNP1 - EPSN1
         else
c
c        le joint se deforme toujours elastiquement en fermeture
c
            DEFP(2) = 0.0D0
c
         endif
C
C     CALCUL DES VARIABLES INTERNES : DEPEQ  STAT  EPOUN LAM1
C
c
c      calcul du delta epsilon plastique entre
c      le pas actuel et le pas precedent
c
         delep1 = defp(1)
         delepn = defp(2)
c
c      produit contracte pour le calcul de la defo. plasti. cumulee
c
         pro1 = ( delep1 * delep1 ) + ( delepn * delepn )
         depeq = sqrt( (2.0D0/3.0D0) * pro1 )
         VARF(1) = EPSEQ + DEPEQ
         VARF(4) = XLAM1 + DELAM1
c
c      calcul de la deformation normale due a l'ouverture seule
c
         if ( abs(snnp11).le.presp ) then
            epoun = delepn
            VARF(2) = EPOUN
         else
            epoun = 0.0d0
            VARF(2) = EPOUN
         endif
c
c     etat du joint
c                stat = 0  levres en contact, avec ou sans cisaillement
c                stat = 1  ouverture du joint pente EC
c                stat = 2  compression du joint pente EC
c                stat = 3  compression du joint pente EF
c
         if ( abs(snnp11).le.presp ) then
            stat = 1.D0
            if ( abs(epnnp1).le.1.D-8 ) then
               stat = 0.D0
            endif
            VARF(3) = STAT
         endif
         if ( (snnp11.gt.sigsei).and.(snnp11.lt.presm) ) then
            stat = 2.D0
            VARF(3) = STAT
         endif
         if (snnp11.le.sigsei) then
            stat = 3.D0
            VARF(3) = STAT
         endif
c
         RETURN
      ELSE
C
         S1NP10=S1NP11
         SNNP10=SNNP11
         S2 = 0.0D0
         CRICIS = CALG(S1NP10,S2,SNNP10,PHI,C)
         CRINOR = SNNP10 - RTRAC
         ITER = ITER + 1
C
         GO TO 2
      ENDIF
C
C
C____________________________________________________________________
C
C
C----------CRITERE DE L'EFFORT NORMAL VIOLE MAIS
C----------CRITERE DE CISAILLEMENT VERIFIE
C
 30   CONTINUE
C
C            ECOULEMENT SELON LE CRITERE F A L'ITERATION ITER
C
      ITER = 1
      DELAM1=0.D0
C
 3    CONTINUE
C
      IF (ITER.GT.MAXITE) THEN
         WRITE(6,*) ' LE CRITERE NORMAL EST VIOLE '
         WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES'
         KERRE = 2
         RETURN
      ENDIF
C
      IF (IMPRI.EQ.1) THEN
         write(6,*) 'PASSAGE N. 3'
         write(6,*) 'CRITERE DE L EFFORT NORMAL VIOLE'
         write(6,*) 'ITERATION',ITER
      ENDIF
C
C
C            CALCUL DU CORRECTEUR PLASTIQUE
C
C
C                DERIVEE DU CRITERE F SUIVANT SIGMA
C
      DFIDS1 = 0.0D0
      DFIDSN = 1.0D0
C
C                MATRICE DE HOOK * DGI/DSIGMA
C
      HODF1 = XMAT(1) * DFIDS1
      HODFN = XMAT(2) * DFIDSN
C
C                FIFI = DFI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
C
      FIFI = (DFIDS1 * HODF1) + (DFIDSN * HODFN)
C
C                RESOLUTION DE L'EQUATION POUR TROUVER LE CORRECTEUR
C                PLASTIQUE DLAM1
C
      IF (ABS(FIFI).LE.PRECIS) THEN
         WRITE(6,*) ' SYSTEME INDETERMINE (F) : ITERATION ',ITER
         KERRE = 2
         RETURN
      ENDIF
C
      DLAM1 = CRINOR / FIFI
      DELAM1 = DELAM1 + DLAM1
      RTRAC =0.D0
C
C                CALCUL DU NOUVEL ETAT DE CONTRAINTES
C
      S1NP11=S1NP10 - (XMAT(1) *  (DLAM1*DFIDS1) )
      SNNP11=SNNP10 - (XMAT(2) *  (DLAM1*DFIDSN) )
C
C   LE TEST DE LA CONVERGENCE
C
      S2 = 0.0D0
      CRICIS = CALG(S1NP11,S2,SNNP11,PHI,C)
      CRINOR = SNNP11 - RTRAC
C
      IF (ABS(CRINOR).LE.PRECIS) THEN
C
C----------LES DEUX CRITERES SONT VERIFIES
C
C
         IF (IMPRI.EQ.1) THEN
            write(6,*) 'LES DEUX CRITERES SONT VERIFIES'
            write(6,*) 'CONVERGENCE A L ITERATION',ITER
         ENDIF
C
C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
C
C
         SIGF(1) = S1NP11
         SIGF(2) = SNNP11
C
         DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1
c
         if ( abs(sigf(2)) .le. presp ) then
c
c            deformation plastique du pas (n+1) moins
c            deformation plastique du pas (n)
c
            DEFP(2) = EPNNP1 - EPSN1
         else
c
c        le joint se deforme toujours elastiquement en fermeture
c
            DEFP(2) = 0.0D0
c
         endif
C
C     CALCUL DES VARIABLES INTERNES : DEPEQ  STAT  EPOUN XLAM1
C
c
c      calcul du delta epsilon plastique entre
c      le pas actuel et le pas precedent
c
         delep1 = defp(1)
         delepn = defp(2)
c
c      produit contracte pour le calcul de la defo. plasti. cumulee
c
         pro1 = ( delep1 * delep1 ) + ( delepn * delepn )
         depeq = sqrt( (2.0D0/3.0D0) * pro1 )
         VARF(1) = EPSEQ + DEPEQ
         VARF(4) = XLAM1 + DELAM1
c
c      calcul de la deformation normale due a l'ouverture seule
c
         if ( abs(snnp11).le.presp ) then
            epoun = delepn
            VARF(2) = EPOUN
         else
            epoun = 0.0d0
            VARF(2) = EPOUN
         endif
c
c     etat du joint
c                stat = 0  levres en contact, avec ou sans cisaillement
c                stat = 1  ouverture du joint pente EC
c                stat = 2  compression du joint pente EC
c                stat = 3  compression du joint pente EF
c
         if ( abs(snnp11).le.presp ) then
            stat = 1.D0
            if ( abs(epnnp1).le.1.D-8 ) then
               stat = 0.D0
            endif
            VARF(3) = STAT
         endif
         if ( (snnp11.gt.sigsei).and.(snnp11.lt.presm) ) then
            stat = 2.D0
            VARF(3) = STAT
         endif
         if (snnp11.le.sigsei) then
            stat = 3.D0
            VARF(3) = STAT
         endif
c
         RETURN
      ELSE
C
         S1NP10=S1NP11
         SNNP10=SNNP11
         S2 = 0.0D0
         CRICIS = CALG(S1NP10,S2,SNNP10,PHI,C)
         CRINOR = SNNP10 - RTRAC
         ITER = ITER + 1
C
         GO TO 3
      ENDIF
C
C
C____________________________________________________________________
C
C
C----------CRITERES DE CISAILLEMENT VIOLE MAIS
C----------CRITERE DE L'EFFORT NORMAL VERIFIE
C
 40   CONTINUE
C
C            ECOULEMENT SELON LE CRITERE G A L'ITERATION ITER
C
      ITER = 1
C
 4    CONTINUE
C
      IF (ITER.GT.MAXITE) THEN
         WRITE(6,*) ' LE CRITERE DE CISAILLEMENT EST VIOLE '
         WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES'
         KERRE = 2
         RETURN
      ENDIF
C
      IF (IMPRI.EQ.1) THEN
         write(6,*) 'PASSAGE N. 4'
         write(6,*) 'CRITERE DE CISAILLEMENT VIOLE'
         write(6,*) 'ITERATION',ITER
      ENDIF
C
C
C            CALCUL DU CORRECTEUR PLASTIQUE
C
C
C                DERIVEE DU CRITERE G SUIVANT SIGMA
C
      if ( s1np10.lt.precis ) then
         dgids1 = -1.d0
      else
         dgids1 =  1.d0
      endif
c
      DZGIDS = TAN(ZMU)
      DGIDSN = TAN(PHI)
C
C                MATRICE DE HOOK * DGI/DSIGMA
C
      HODG1 = XMAT(1) * DGIDS1
      HODGN = XMAT(2) * DZGIDS
C
C                GIGI = DGI/DSIGMA * MATRICE DE HOOK * DGI/DSIGMA
C
      GIGI = (DGIDS1 * HODG1) + (DGIDSN * HODGN)
C
C                RESOLUTION DE L'EQUATION POUR TROUVER LE CORRECTEUR
C                PLASTIQUE DLAM2
C
      IF (ABS(GIGI).LE.PRECIS) THEN
         WRITE(6,*) ' SYSTEME INDETERMINE (G) : ITERATION ',ITER
         KERRE = 2
         RETURN
      ENDIF
C
      DLAM2 = CRICIS / GIGI
C
C                CALCUL DU NOUVEL ETAT DE CONTRAINTES
C
      S1NP11=S1NP10 - (XMAT(1) *  (DLAM2*DGIDS1) )
      SNNP11=SNNP10 - (XMAT(2) *  (DLAM2*DZGIDS) )
c
      if (SNNP11.GE.0.0d0) THEN
         SNNP11 = 0.0d0
      ENDIF
cccssszzzeeerrr       if (abs(zmu).le.presp) then
cccssszzzeeerrr          differ = abs( depsn ) - epsn1
cccssszzzeeerrr          differ = abs(differ)/abs( depsn )
cccssszzzeeerrr          if (depsn.lt.preem) then
cccssszzzeeerrr             if (epnnp1.gt.preem) then
cccssszzzeeerrr                SNNP11 = 0.0d0
cccssszzzeeerrr             endif
cccssszzzeeerrr          else
cccssszzzeeerrrc
cccssszzzeeerrr             SNNP11 = 0.0d0
cccssszzzeeerrr          endif
cccssszzzeeerrr       endif
C
C      TEST DE LA CONVERGENCE
C
      S2 = 0.0D0
      CRICIS = CALG(S1NP11,S2,SNNP11,PHI,C)
      CRINOR = SNNP11 - RTRAC
C
      IF (ABS(CRICIS).LE.PRECIS) THEN
C
         IF (INDIC.EQ.1) THEN
C
C----------1ER TRAJET : ON PARCOURT ALP * DEPSN
C
            DEFPN = -epin0(2)
            S1INT = S1NP11
            DEFP11 = EP1NP1 - ( DD(1,1) * S1NP11 )
            GO TO 98
         ELSE
            IF (IFIN.EQ.0) THEN
C
C----------LES DEUX CRITERES SONT VERIFIES
C
C
               IF (IMPRI.EQ.1) THEN
                  write(6,*) 'LES DEUX CRITERES SONT VERIFIES'
                  write(6,*) 'CONVERGENCE A L ITERATION',ITER
               ENDIF
C
C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
C
C
               SIGF(1) = S1NP11
               SIGF(2) = SNNP11
C
               DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1
c
               if ( abs(sigf(2)) .le. presp ) then
c
c            deformation plastique du pas (n+1) moins
c            deformation plastique du pas (n)
c
                  DEFP(2) = EPNNP1 - EPSN1
               else
c
c        le joint se deforme toujours elastiquement en fermeture
c
                  DEFP(2) = 0.0D0
c
               endif
            ELSE
C
C----------CAS DU PASSAGE DIRECT D'UN ETAT OUVERT A UN ETAT COMPRIME
C
               DEFP12  = EP1NP1 - (DD(1,1)*S1NP11)
               DEFP(1) =  DEFP12 - epin0(1)
               DEFP(2) = -epin0(2)
               SIGF(1) = S1NP11
               SIGF(2) = SNNP11
               IFIN    = 0
            ENDIF
C
         ENDIF
C
C     CALCUL DES VARIABLES INTERNES : DEPEQ  STAT  EPOUN XLAM1
C
c
c      calcul du delta epsilon plastique entre
c      le pas actuel et le pas precedent
c
         delep1 = defp(1)
         delepn = defp(2)
c
c      produit contracte pour le calcul de la defo. plasti. cumulee
c
         pro1 = ( delep1 * delep1 ) + ( delepn * delepn )
         depeq = sqrt( (2.0D0/3.0D0) * pro1 )
         VARF(1) = EPSEQ + DEPEQ
         VARF(4) = XLAM1
c
c      calcul de la deformation normale due a l'ouverture seule
c
         if ( abs(snnp11).le.presp ) then
            epoun = delepn
            VARF(2) = EPOUN
         else
            epoun = 0.0d0
            VARF(2) = EPOUN
         endif
c
c     etat du joint
c                stat = 0  levres en contact, avec ou sans cisaillement
c                stat = 1  ouverture du joint pente EC
c                stat = 2  compression du joint pente EC
c                stat = 3  compression du joint pente EF
c
         if ( abs(snnp11).le.presp ) then
            stat = 1.D0
            if ( abs(epnnp1).le.1.D-8 ) then
               stat = 0.D0
            endif
            VARF(3) = STAT
         endif
         if ( (snnp11.gt.sigsei).and.(snnp11.lt.presm) ) then
            stat = 2.D0
            VARF(3) = STAT
         endif
         if (snnp11.le.sigsei) then
            stat = 3.D0
            VARF(3) = STAT
         endif
c
         RETURN
      ELSE
C
         S1NP10=S1NP11
         SNNP10=SNNP11
         S2 = 0.0D0
         CRICIS = CALG(S1NP10,S2,SNNP10,PHI,C)
         CRINOR = SNNP10 - RTRAC
         ITER = ITER + 1
C
         GO TO 4
      ENDIF
C
C
C____________________________________________________________________
C
      RETURN
      END