C RIC3NL    SOURCE    OF166741  25/11/04    21:16:04     12349          
C sub ricnl3d
       SUBROUTINE RIC3NL(wrk52,wrk53,wrk54,nvari,iecou)
C
C====&===1=========2=========3=========4=========5=========6=========7==
C Commentaires : Subroutine permettant de mettre en oeuvre le
C                    modele RICRAG (F. Ragueneau & B. Richard) pour représenter
C                le comportement 3D d'un béton ordinaire sous
C                sollicitations complexes en non local
C Traits       : - Endommagement scalaire
C                - Glissement frottant
C                - 3D
C                - Monotone et cyclique
C                - Non localité portant sur le seuil d'endommagement (Fd)
C                - Bon comportement unilatteral
C Auteurs      : B. Richard (doctorant) & F. Ragueneau (LCPC, FDOA & BCC)
C Date         : Mai 2008
C====&===1=========2=========3=========4=========5=========6=========7==
C
C-----DECLARATION GENERALE----------------------------------------------
C
      IMPLICIT REAL*8(A-H,O-Z)
C
C-----LISTE DES INCLUDES A CHARGER--------------------------------------
C

-INC PPARAM
-INC CCOPTIO
-INC DECHE
-INC CCREEL
C
C-----DECLARATION GENERALE----------------------------------------------
C
C       REAL*8 XMAT(9),VAR0(30),VARF(30),SIGF(6),DEPST(6)
      REAL*8 XZER,UN,DEUX,EPSI(3,3)
      REAL*8 E,NU,FT,ALIND,GAMMA1,A1,ALDIR
      REAL*8 EPS0,K,G,Y0,DEPSDEV(3,3),DTRA
      REAL*8 UNIT(3,3),DEPSI(3,3),EPSTFBIS(6)
      REAL*8 EPSIPI(3,3),ALPHA(3,3),Z,D,DONMDNL,DONMINL,FFFZ,DOLD,FORCED
      REAL*8 FORCEI,NRUPT,SEUILT,DEVPI(3,3)
      REAL*8 TRA,EPSDEV(3,3),SIGMA(3,3),VPEPS(3),SEUIL1NL
      REAL*8 YD,YN,EPSDIR(3,3),EPSIND(3,3),Dn
      REAL*8 SIGMPI(3,3),SIGMOD(3,3),SIGMOX(3,3),X(3,3),TERME6(3,3)
      REAL*8 TERME7(3,3),TERME8(3,3),TERME9(3,3),TROGO,DEVTRO(3,3)
      REAL*8 DEVDGS(3,3)
C       CHARACTER*8 CMATE
C
C-----DEFINITION DE PARAMETRES------------------------------------------
C
      PARAMETER (XZER=0.0D0 , UN=1.0D0 , DEUX=2.0D0, XCOEF=0.00D0)
C
C-----OUVERTURE DE SEGMENTS DE TRAVAIL----------------------------------
C
C-----MISE EN DONNEE----------------------------------------------------
C
      CMATE='ISOTROPE'
      E      = XMAT(1)
      NU     = XMAT(2)
      FT     = XMAT(5)
      ALIND  = XMAT(6)
      GAMMA1 = XMAT(7)
      A1     = XMAT(8)
      ALDIR  = XMAT(9)

      EPS0   = FT/E
      K      = E/(1.D0-2.D0*NU)
      G      = E/(2.D0*(1.D0+NU))
      Y0     = (K/6.0D0*EPS0**2.D0)

        DO J=1,3
         DO I=1,3
            IF (I.EQ.J) THEN
               UNIT(I,J)=UN
            ELSE
               UNIT(I,J)=XZER
            ENDIF
         ENDDO
      ENDDO
C
C-----PASSAGE DES VARIABLES INTERNES------------------------------------
C
C     GLISSEMENT (EPSILON_PI)
C
      EPSIPI(1,1)=VAR0(2)
      EPSIPI(2,2)=VAR0(3)
      EPSIPI(3,3)=VAR0(4)

      EPSIPI(1,2)=VAR0(5)
      EPSIPI(1,3)=VAR0(6)
      EPSIPI(2,3)=VAR0(7)

      EPSIPI(2,1)=VAR0(5)
      EPSIPI(3,1)=VAR0(6)
      EPSIPI(3,2)=VAR0(7)
C
C     ECROUISSAGE CINEMATIQUE (ALPHA)
C
      ALPHA(1,1)=VAR0(8)
      ALPHA(2,2)=VAR0(9)
      ALPHA(3,3)=VAR0(10)
C
      ALPHA(1,2)=VAR0(11)
      ALPHA(1,3)=VAR0(12)
      ALPHA(2,3)=VAR0(13)
C
      ALPHA(2,1)=VAR0(11)
      ALPHA(3,1)=VAR0(12)
      ALPHA(3,2)=VAR0(13)
C
C     ENDOMMAGEMENT (D)
C
      D=VAR0(14)
C
C     FORCE ECROUISSAGE ISOTROPE
C
      FFFZ=VAR0(17)
C
C     DEFORMATION TOTALE MISE DANS EPSI
C
      DO I=1,6
           EPSTFBIS(I)=VAR0(18+I)+DEPST(I)
      ENDDO

      DO I=1,3
           EPSI (I,I)=VAR0(18+I)+DEPST(I)
           DEPSI(I,I)=DEPST(I)
      ENDDO

      EPSI(1,2)=VAR0(22)+DEPST(4)*0.5D0
      EPSI(1,3)=VAR0(23)+DEPST(5)*0.5D0
      EPSI(2,3)=VAR0(24)+DEPST(6)*0.5D0
      EPSI(2,1)=VAR0(22)+DEPST(4)*0.5D0
      EPSI(3,1)=VAR0(23)+DEPST(5)*0.5D0
      EPSI(3,2)=VAR0(24)+DEPST(6)*0.5D0

      DEPSI(1,2)=DEPST(4)*0.5D0
      DEPSI(1,3)=DEPST(5)*0.5D0
      DEPSI(2,3)=DEPST(6)*0.5D0
      DEPSI(2,1)=DEPST(4)*0.5D0
      DEPSI(3,1)=DEPST(5)*0.5D0
      DEPSI(3,2)=DEPST(6)*0.5D0

      SIGMA(1,1)=SIG0(1)
      SIGMA(2,2)=SIG0(2)
      SIGMA(3,3)=SIG0(3)
      SIGMA(1,2)=SIG0(4)
      SIGMA(1,3)=SIG0(5)
      SIGMA(2,3)=SIG0(6)
      SIGMA(2,1)=SIG0(4)
      SIGMA(3,1)=SIG0(5)
      SIGMA(3,2)=SIG0(6)
C
C     ECROUISSAGE
C
      FORCED = VAR0(25)
      FORCEI = VAR0(26)
C
C     MATRICE DE HOOKE
C
      AUX0=E/((UN+NU)*(UN-NU-NU))
      AUX=AUX0*(UN-NU)
      AUX1=AUX0*NU
      XGEGE=0.5D0*E/(UN+NU)

      DDHOOK(1,1)=AUX
      DDHOOK(2,1)=AUX1
      DDHOOK(3,1)=AUX1

      DDHOOK(1,2)=AUX1
      DDHOOK(2,2)=AUX
      DDHOOK(3,2)=AUX1

      DDHOOK(1,3)=AUX1
      DDHOOK(2,3)=AUX1
      DDHOOK(3,3)=AUX

      DDHOOK(4,4)=XGEGE
      DDHOOK(5,5)=XGEGE
      DDHOOK(6,6)=XGEGE
C-----------------------------------------------------------------------
C-----DEBUT INTEGRATION NUMERIQUE---------------------------------------
C-----------------------------------------------------------------------
      TRA     = (EPSI(1,1)+EPSI(2,2)+EPSI(3,3))/3.D0
      DTRA    = (DEPSI(1,1)+DEPSI(2,2)+DEPSI(3,3))/3.D0

        DO J=1,3
         DO I=1,3
            EPSDEV(I,J)  = EPSI(I,J)-TRA*UNIT(I,J)
            DEPSDEV(I,J) = DEPSI(I,J)-DTRA*UNIT(I,J)
         ENDDO
      ENDDO

        DO J=1,3
         DO I=1,3
                SIGMA(I,J)=SIGMA(I,J)+K*DTRA+2.0D0*G*DEPSDEV(I,J)
         ENDDO
      ENDDO
C
C-----ENDOMMAGEMENT------------------------------------------------------
C
C       CALCUL DES VALEURS PRINCIPALES
C
        CALL JACOD3(EPSI ,3,VPEPS)
C
C       CALCUL DE EPSILON DIRECT
C
        IF (TRA.GE.XZER) THEN
         XPDC=1.0D0
c          GAMMA1 = 7.0D8
c          A1     = 7.0D-5
      ELSE
         XPDC=XZER
      ENDIF

        DO J=1,3
         DO I=1,3
            IF (I.EQ.J) THEN
               EPSDIR(I,J)=0.5D0*(ABS(VPEPS(I))+VPEPS(I))*XPDC
            ELSE
               EPSDIR(I,J)=XZER
            ENDIF
         ENDDO
      ENDDO
C
C       CALCUL DE EPSILON INDUIT
C
        DO J=1,3
         DO I=1,3
            IF (I.EQ.J) THEN
               EPSIND(I,J)=0.5D0*(ABS(VPEPS(I))+VPEPS(I))-EPSDIR(I,J)
            ELSE
               EPSIND(I,J)=XZER
            ENDIF
         ENDDO
      ENDDO
C
C       TAUX D ENERGIES CORRESPONDANTES
C
      DONMDNL=XZER
      DONMINL=XZER

      DO J=1,3
           DO I=1,3
              DONMDNL=DONMDNL+EPSDIR(I,J)*EPSDIR(I,J)
              DONMINL=DONMINL+EPSIND(I,J)*EPSIND(I,J)
           ENDDO
      ENDDO

        SEUIL1NL=(K/6.0D0) * (DONMDNL*XPDC+DONMINL*(1.0D0-XPDC))

      IF (ISTEP.EQ.0) THEN
         SEUIL1=SEUIL1NL
         VARF(1)=SEUIL1
      ELSE IF (ISTEP.EQ.1) THEN
         VARF(1)=SEUIL1NL
         DO I=2,26
          VARF(I)=VAR0(I)
         ENDDO
         GOTO 2000
      ELSE IF (ISTEP.EQ.2) THEN
         SEUIL1=VAR0(1)
         VARF(1)=SEUIL1
      ENDIF
C
C       CALCUL NON LOCAL - ON MOYENNE LE SEUIL
C
      FFFZ=FORCED*XPDC+FORCEI*(1.0D0-XPDC)
        SEUILT=SEUIL1-(FFFZ+Y0)
C
C       EVOLUTION ENDOMMAGEMENT
C
        IF (SEUILT.GT.XZER) THEN
           IF (TRA.GE.XZER) THEN
              Dn      = 1.0D0-1.0D0/(1.0D0+ALDIR*(SEUIL1-Y0))
              FORCED = SEUIL1-Y0
           ELSE
              Dn      = 1.0D0-1.0D0/(1.0D0+ALIND*(SEUIL1-Y0))
              FORCEI = SEUIL1-Y0
           ENDIF
C
C     ON VERIFIE D COISSANT
C
        IF (Dn.GE.D) THEN
           D=Dn
        ENDIF
C
C     ON LIMITE LA VALEUR DE D
C
      IF (D.GE.0.99999D0) THEN
         D=0.99999D0
      ENDIF
C
C       ON BLOQUE LES ECROUISSAGES
C
           FFFZ=FORCED*XPDC+FORCEI*(1.0D0-XPDC)
        ENDIF

      DO J=1,3
           DO I=1,3
              X(I,J)=GAMMA1*ALPHA(I,J)
           ENDDO
      ENDDO
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C                               SUBSTEPPING
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        NENTIER=100
        DO KIN=1,NENTIER

        DO I=1,6
           EPSTFBIS(I)=VAR0(18+I)+KIN*DEPST(I)/NENTIER
      ENDDO

        EPSI(1,1)=EPSTFBIS(1)
        EPSI(2,2)=EPSTFBIS(2)
        EPSI(3,3)=EPSTFBIS(3)
        EPSI(1,2)=EPSTFBIS(4)/2.0D0
        EPSI(1,3)=EPSTFBIS(5)/2.0D0
        EPSI(2,3)=EPSTFBIS(6)/2.0D0
        EPSI(2,1)=EPSTFBIS(4)/2.0D0
        EPSI(3,1)=EPSTFBIS(5)/2.0D0
        EPSI(3,2)=EPSTFBIS(6)/2.0D0

        TRA   = (EPSI(1,1)+EPSI(2,2)+EPSI(3,3))/3.D0
        TRAPI = (EPSIPI(1,1)+EPSIPI(2,2)+EPSIPI(3,3))/3.D0

        DO J=1,3
          DO I=1,3
            EPSDEV(I,J)=EPSI(I,J)-TRA*UNIT(I,J)
             DEVPI(I,J)=EPSIPI(I,J)-TRAPI*UNIT(I,J)
          ENDDO
      ENDDO

        DO J=1,3
         DO I=1,3
                SIGMA(I,J)  = (1.0D0-D)*K*TRA*UNIT(I,J)+
     &                    (1.0D0-D)*2.0D0*G*EPSDEV(I,J)+
     &                    K*D*(TRA-TRAPI)*UNIT(I,J)+
     &                2.0D0*G*D*(EPSDEV(I,J)-DEVPI(I,J))
         ENDDO
      ENDDO

         DO J=1,3
          DO I=1,3
            SIGMPI(I,J)=K*D*(TRA-TRAPI)*UNIT(I,J)+
     &              2.0D0*G*D*(EPSDEV(I,J)-DEVPI(I,J))
          ENDDO
       ENDDO
C
C====&===1=========2=========3=========4=========5=========6=========7==
C-----DEBUT INTEGRATION GLISSEMENT--------------------------------------
C====&===1=========2=========3=========4=========5=========6=========7==
C
        DO 500 IREMP=1,500
C
C     1) EVALUATION DU SEUIL
C     ----------------------
        DO 101 J=1,3
          DO 201 I=1,3
            SIGMOX(I,J)=XZER
            SIGMOX(I,J)=SIGMPI(I,J)-X(I,J)
 201      CONTINUE
 101    CONTINUE


        TRASX= (SIGMOX(1,1)+SIGMOX(2,2)+SIGMOX(3,3))/3.D0

        DO 102 J=1,3
          DO 202 I=1,3
            SIGMOD(I,J)=XZER
            SIGMOD(I,J)=SIGMOX(I,J)-TRASX*UNIT(I,J)
 202      CONTINUE
 102    CONTINUE

        TERME4=XZER
        DO 103 J=1,3
          DO 203 I=1,3
            TERME4=TERME4+SIGMOD(I,J)*SIGMOD(I,J)
 203      CONTINUE
 103    CONTINUE
C
C         CALCUL DU J2(SIGMPI-X)
C
        TRASIN = SIGMA(1,1)+SIGMA(2,2)+SIGMA(3,3)
        TRASIN = 0.5D0*(-ABS(TRASIN)+TRASIN)*0.0D0

        TERME4=SQRT(1.5D0*TERME4)+(TRASIN/3.0D0)*0.80D0

        SEUIL2=TERME4+XCOEF*TRASX
        IF (IREMP.EQ.1) THEN
          SEUINI=SEUIL2
          DLAM2=XZER
          COCO=2
        ENDIF

        IF (SEUINI.GT.1.0D0) THEN
          CRIT=SEUIL2/SEUINI
        ELSE
          CRIT=XZER
        ENDIF
C
C       2) ON VERIFIE LA VALEUR RELATIVE DU SEUIL
C       -----------------------------------------
        IF ((CRIT.LE.1.0D-5).OR.(DLAM2.LE.1.0D-10.AND.IREMP.GT.1).OR.
     &       (D.EQ.XZER))THEN
          GOTO 666
        ELSE
        IF (IREMP.EQ.1) COCO=1
C
C       3) CALCUL DES DERIVÉES POUR LA REDESCENTE
C       -----------------------------------------
        DO J=1,3
          DO I=1,3
            TERME6(I,J)=XZER
            TERME7(I,J)=XZER
            TERME8(I,J)=XZER
            TERME9(I,J)=XZER
            ENDDO
          ENDDO
C
C       CALCUL DE DF/DS....................................................
C
        DO J=1,3
          DO I=1,3
            TERME6(I,J)=1.5D0*SIGMOX(I,J)/TERME4+XCOEF*1.0D0/3.D0*
     &                  UNIT(I,J)
            ENDDO
          ENDDO
C
C       CALCUL DE DF/DX....................................................
C
        DO J=1,3
          DO I=1,3
            TERME7(I,J)=-TERME6(I,J)
            ENDDO
          ENDDO
C
C       CALCUL DE DG/DS....................................................
C
        DO J=1,3
          DO I=1,3
            TERME8(I,J)=1.5D0*SIGMOX(I,J)/TERME4+XCOEF*1.0D0/3.D0*
     &                  UNIT(I,J)
            ENDDO
          ENDDO
C
C       CALCUL DE DG/DX....................................................
C
        DO J=1,3
          DO I=1,3
            TERME9(I,J)=-1.5D0*SIGMOX(I,J)/TERME4+A1*X(I,J)-
     &                  XCOEF*1.0D0/3.D0*UNIT(I,J)
            ENDDO
          ENDDO
C
C       4) CALCUL DU MULTIPLICATEUR DE GLISSEMENT
C       -----------------------------------------
          TRDGZ=(TERME8(1,1)+TERME8(2,2)+TERME8(3,3))/3.D0

          DO J=1,3
          DO I=1,3
            DEVDGS(I,J)=TERME8(I,J)-TRDGZ*UNIT(I,J)
            ENDDO
          ENDDO

        TERM10=XZER
        DO J=1,3
          DO I=1,3
            TERM10=TERM10+TERME6(I,J)*(K*TRDGZ*UNIT(I,J)+
     &             2.0D0*G*DEVDGS(I,J))
            ENDDO
          ENDDO

        TERM11=XZER
        DO J=1,3
          DO I=1,3
            TERM11=TERM11+TERME7(I,J)*TERME9(I,J)
            ENDDO
          ENDDO

        DLAM2=SEUIL2/(TERM10+GAMMA1*TERM11)
C
C       5) ACTUALISATION DES VARIABLES FORCES
C       -------------------------------------
        DO J=1,3
          DO I=1,3
            SIGMPI(I,J)=SIGMPI(I,J)-DLAM2*D*(K*TRDGZ*UNIT(I,J)+
     &             2.0D0*G*DEVDGS(I,J))
            ENDDO
          ENDDO

        DO J=1,3
          DO I=1,3
            X(I,J)=X(I,J)-GAMMA1*DLAM2*TERME9(I,J)
            ENDDO
          ENDDO

        ENDIF
 500  CONTINUE
C
C-----ACTUALISATION DES VARIABLES FLUX (S'IL Y A LIEU)------------------
C
C     1) ALPHA :
C       -------
 666  CONTINUE
      IF (COCO.EQ.1) THEN
        DO 303 J=1,3
          DO 403 I=1,3
            ALPHA(I,J)=X(I,J)/GAMMA1
 403      CONTINUE
 303    CONTINUE
C
C     2) EPSILON_PI :
C        ------------
        dplus=max(d,xpetit)
        TRA_Z = (SIGMPI(1,1)+SIGMPI(2,2)+SIGMPI(3,3))/3.D0
        TROGO = TRA - TRA_Z/(Dplus*K)
        DO J=1,3
          DO I=1,3
            DEVTRO(I,J)=EPSDEV(I,J)-(SIGMPI(I,J)-TRA_Z*UNIT(I,J))
     &                                /(2.0D0*G*Dplus)
            ENDDO
        ENDDO

        DO J=1,3
          DO I=1,3
            EPSIPI(I,J)=TROGO*UNIT(I,J)+DEVTRO(I,J)
            ENDDO
        ENDDO

      ENDIF
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        ENDDO
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C-----CORRECTION DES CONTRAINTES
C
      TRAPOS = ( ABS(TRA)+TRA)*0.5D0
      TRANEG = (-ABS(TRA)+TRA)*0.5D0

       DO J=1,3
         DO I=1,3
              SIGMA(I,J)  = (1.0D0-D)*K*TRAPOS*UNIT(I,J)+
     &                              K*TRANEG*UNIT(I,J)+
     &                          (1.0D0-D)*2.0D0*G*EPSDEV(I,J)+
     &                           SIGMPI(I,J)
         ENDDO
      ENDDO
C------------------------------------------------------------------------
C-----FIN INTEGRATION NUMERIQUE------------------------------------------
C------------------------------------------------------------------------
      SIGF(1)= SIGMA(1,1)
      SIGF(2)= SIGMA(2,2)
      SIGF(3)= SIGMA(3,3)
      SIGF(4)= (SIGMA(1,2)+SIGMA(2,1))*0.5D0
      SIGF(5)= (SIGMA(1,3)+SIGMA(3,1))*0.5D0
      SIGF(6)= (SIGMA(2,3)+SIGMA(3,2))*0.5D0
C
C     GLISSEMENT (EPSILON_PI)
C
      VARF(2)=EPSIPI(1,1)
      VARF(3)=EPSIPI(2,2)
      VARF(4)=EPSIPI(3,3)
      VARF(5)=EPSIPI(1,2)
      VARF(6)=EPSIPI(1,3)
      VARF(7)=EPSIPI(2,3)
C
C     ECROUISSAGE CINEMATIQUE (ALPHA)
C
      VARF(8 )=ALPHA(1,1)
      VARF(9 )=ALPHA(2,2)
      VARF(10)=ALPHA(3,3)
      VARF(11)=ALPHA(1,2)
      VARF(12)=ALPHA(1,3)
      VARF(13)=ALPHA(2,3)
C
C     ENDOMMAGEMENT (D)
C
      VARF(14)=D
C
C     FORCE LIEE A L ECROUISSAGE ISOTROPE
C
      VARF(17)=FFFZ
C
C     DEFORMATION TOTALE
C
      DO I=1,6
         VARF(18+I)=EPSTFBIS(I)
        ENDDO
C
C     ECROUISSAGES
C
        VARF(25)=FORCED
        VARF(26)=FORCEI
C
C     CONTRAINTES DE FROTTEMENT
C
*      VARF(27)=SIGMPI(1,1)
*      VARF(28)=SIGMPI(2,2)
*      VARF(29)=SIGMPI(3,3)
*      VARF(30)=SIGMPI(1,3)

 2000 CONTINUE

        RETURN
        END