Sources Esope | Cast3M

cendom
C CENDOM    SOURCE    OF166741  25/11/04    21:15:13     12349          
       SUBROUTINE CENDOM(wrk52,wrk53,wrk54,WRK6,WRK7,WRK8,WRK9,
     1           NVARI, T0,T1,TREF,IB,IGAU,IFOUR7,NCOUR7,iecou)
 
c NSTRSS &lt;= NSTRS , MFR1 &lt;- MFR , XCARB &lt;- XCAR
C ==================================================================
C CE SOUS-PROGRAMME EST APPELE DANS "PLAST".
C C'EST UN INTEGRATEUR GENERAL DES LOIS DE COMPORTEMENT DES
C MATERIAUX ENDOMMAGEABLES EN REGIME ANISOTHERME.
C
C ENTREES:
C -------
C WRK0,WR00,WRK1,WRK6,WRK7,WRK8,WRK9,WTRAV = SEGMENTS DE TRAVAIL
C NSTRSS          = NBR. DE COMPOSANTES DES CONTR. OU DES DEFORM.
C SIG0(NSTRSS)    = CONTR. AU DEBUT DU PAS D'INTEGRATION
C DEPST(NSTRSS)   = INCREMENT DES DEFORM. ELASTIQUES
C NVARI          = NBR. DE VARIABLES INTERNES
C VAR0(NVARI)    = VARIABLES INTERNES AU DEBUT DU PAS D'INTEGRATION
C                  CE TABLEAU CONTIENT DANS L'ORDRE:
C                - LA DEFORM. PLASTIQUE CUMULEE
C                - LES VARIABLES INTERNES PILOTANT LES EQUATIONS DU
C                  MODELE
C                - LES DEFORM. PLASTIQUES
C MFR1            = INDICE DE LA FORMULATION MECANIQUE; SEULEMENT
C                  MASSIF OU COQUE POUR LES MATERIAUX ENDOMMAGEABLES
* CMATE = NOM DU MATERIAU
* N2EL  = NBRE D ELEMENTS DANS SEGMENT DE HOOKE
* N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE
* MFR1   = NUMERO DE LA FORMULATION
* IFOUR7= OPTION DE CALCUL
* IB    = NUMERO DE L ELEMENT COURANT
* IGAU  = NUMERO DU POINT COURANT
* EPAIST= EPAISSEUR
* NBPGAU= NBRE DE POINTS DE GAUSS
* MELE  = NUMERO DE L ELEMENT FINI
* NPINT = NBRE DE POINTS D INTEGRATION
* NBGMAT= NBRE DE POINTS DANS SEGMENT DE CARACTERISTIQUES
* NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES
* SECT  = SECTION
* LHOOK = TAILLE DE LA MATRICE DE HOOKE
C ICARA          = NBR. DE CARACT. GEOMETRIQUES DES ELEMENTS FINIS
C XCARB(ICARA)    = CARACT. GEOMETRIQUES DES ELEMENTS FINIS
C INPLAS         = INDICE DU MODELE DECRIVANT LE MATERIAU ENDOMMAGEABLE
C NYKX          = DIMENSION DE YKX
C YKX(NYKX)   = TABLEAU CONTENANT LES COURBES DE TRACTION(EPS,SIGMA)
C                  A DES TEMPERATURES DONNEES (Ti);IL EST DE LA FORME
C                  SUIVANTE :
C                  (T1,X1,F(X1,T1),X2,F(X2,T1),..,T2,Y1,F(Y1,T2),Y2,..)
C NNKX           = NBRE DE COURBES DE TRACTION
C NKX            = TABLEAU CONTENANT LE NBRE DE PTS POUR CHAQUE COURBE
C                  DE TRACTION (EX:SI NKX(1)=4,LA COURBE DE TRACTION
C                  A LA TEMPERATURE T1 SERA DEFINIE PAR 4 POINTS)
C PRECIS         = PRECISION DES ITERATIONS INTERNES
C TO             = TEMPERATURE INITIALE
C T1             = TEMPERATURE FINALE
C Tref           = TEMPERATURE DE REFERENCE
C NYOG          = 2*NBRE DE POINTS DE LA COURBE T-->YOG(T)
C YOG(NYOG)    = TABLEAU CONTENANT LES VALEURS CONNUES DU MODULE
C                  D'YOUNG E A DES TEMPERATURES PRECISES T
C NYNU           = 2*NBRE DE POINTS DE LA COURBE T-->YNU(T)
C YNU(NYNU)      = TABLEAU CONTENANT LES VALEURS CONNUES DU COEFF. DE
C                  POISSON YNU A DES TEMPERATURES PRECISES T
C NYALFA          = 2*NBRE DE POINTS DE LA COURBE T-->YALFA(T)
C YALFA(NYALFA)    = TABLEAU CONTENANT LES VALEURS CONNUES DU COEFF. DE
C                  DILAT. THERMIQUE YALFA A DES TEMPERATURES PRECISES T
C NYN          = 2*NBRE DE POINTS DE LA COURBE T-->YN(T)
C YN(NYN)    = TABLEAU CONTENANT LES VALEURS CONNUES DU SEUIL
C                  D'ENDOMT YN A DES TEMPERATURES PRECISES T
C NYM            = 2*NBRE DE POINTS DE LA COURBE T-->YM(T)
C YM(NYM)        = TABLEAU CONTENANT LES VALEURS CONNUES DE L'ENDOMT
C                  CRITIQUE YM A DES TEMPERATURES PRECISES T
C NYKK          = 2*NBRE DE POINTS DE LA COURBE T-->YKK(T)
C YKK(NYKK)    = TABLEAU CONTENANT LES VALEURS CONNUES DE LA DEFORM.
C                  A RUPTURE YKK A DES TEMPERATURES PRECISES T
C
C SORTIES:
C -------
C NCOUR7     = DIMENSION DE TRUC
C TRUC       = TABLEAU DE DIMENSION LA PLUS GRANDE PARMI CELLES
C              DEFINIES DANS NKX;IL CONTIENDRA LA COURBE DE
C              TRACTION A LA TEMPERATURE T,OBTENUE GRACE AU SOUS-
C              PROGRAMME TRACTT.
C SIGF(NSTRSS)= CONTR. A LA FIN DU PAS D'INTEGRATION
C VARF(NVARI)= VARIABLES INTERNES A LA FIN DU PAS D'INTEGRATION
C DEFP(NSTRSS)= INCREMENT DES DEFORM. PLASTIQUES A LA FIN DU PAS
C              D'INTEGRATION
C KERRE      = INDICE QUI REGIT LES ERREURS
C            = 77 SI LA DEFORM. PLAST. CUMULEE ENDOMMAGEE (2IEME VAR.
C                 INT.) EST EN DEHORS DE LA COURBE DE TRACTION, DS.
C                 LE CAS DE L'ECROUISSAGE ET DE L'ENDOMM. ISOTROPES.
C                 CECI PEUT SE PRODUIRE SUITE A L'APPEL A "CRIDAM"
C            = 99 SI LA FORMULATION MECANIQUE N'EST PAS DISPONIBLE
C                  POUR LE MODELE CONSIDERE OU S'IL Y A INCOMPATIBILITE
C                 ENTRE MFR1 ET IFOUR7
C ==================================================================
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC DECHE
 
-INC TECOU
 
      SEGMENT WRK6
        REAL*8 BB(NSTRS,NNVARI),R(NSTRS),XMU(NSTRS)
        REAL*8 S(NNVARI),QSI(NNVARI),DDR(NSTRS),BBS(NSTRS)
        REAL*8 SIGMA(NSTRS),SIGGD(NSTRS)
      ENDSEGMENT
 
      SEGMENT WRK7
        REAL*8 F(NCOUR7,2),W(NCOUR7),TRUC(NCOUR7)
      ENDSEGMENT
 
      SEGMENT WRK8
        REAL*8 DD(NSTRS,NSTRS),DDV(NSTRS,NSTRS),DDINV(NSTRS,NSTRS)
        REAL*8 DDINVp(NSTRS1,NSTRS1)
      ENDSEGMENT
 
      SEGMENT WRK9
        REAL*8 YOG(NYOG),YNU(NYNU),YALFA(NYALFA),YSMAX(NYSMAX)
        REAL*8 YN(NYN),YM(NYM),YKK(NYKK),YALFA1(NYALF1)
        REAL*8 YBETA1(NYBET1),YR(NYR),YA(NYA),YKX(NYKX),YRHO(NYRHO)
        REAL*8 SIGY(NSIGY)
        INTEGER NKX(NNKX)
      ENDSEGMENT
 
      DIMENSION DEPS(6)
 
C ------------------------------------------------
C DEFINITION DU NBR. MAXIMUM D'ITERATIONS INTERNES
C ------------------------------------------------
      DATA ITMAX/30/
 
      mfrl = iecou.MFR1
      nstl = iecou.NSTRSS
 
      NXMATT = NMATT
 
      ZEERO=0.D0
      PRECI=PRECIS
      PRECIG=-2.D0*PRECIS
      PRECG=ABS(PRECIG)
      INDIC=0
      IALPH=0
      T=T1
      YUNGV=0.D0
      XNUV=0.D0
C
C ================================================================
C MAJ. POUR LE CAS ELASTIQUE ET TEST POUR SAVOIR SI ON A PLASTIFIE
C ================================================================
      CALL ZDANUL(DEFP,nstl)
      ALPH=1.D0
      ALPH0=0.D0
      ALPH1=1.D0
      X=VAR0(2)
 
      CALL INITTE(NYOG,YOG,NYNU,YNU,NYRHO,YRHO,NYALFA,YALFA,NYN,
     1     YN,NYM,YM,NYKK,YKK,NYKX,YKX,NNKX,NKX,X,T0,
     2                  XMAT1,NXMATT,TRUC,NCOUR7)
      CALL INITTE(NYOG,YOG,NYNU,YNU,NYRHO,YRHO,NYALFA,YALFA,NYN,
     1     YN,NYM,YM,NYKK,YKK,NYKX,YKX,NNKX,NKX,X,T1,
     2                  XMAT2,NXMATT,TRUC,NCOUR7)
      ALFA0=XMAT1(4)
      ALFA1=XMAT2(4)
      CTEPS=ALFA0*(T0-TREF)-ALFA1*(T1-TREF)
 
C==================================================
C  CALCUL DE L INCREMENT DE CONTRAINTES
C  DEPST EST L INCREMENT DE DEFORMATIONS ELASTIQUES
C==================================================
C
      CALL ELAST1(1,IFOUR,VAR0,NNVARI,XMAT2,NXMATT,YUNGV,XNUV,
     1          XCARB,ICARA,mfrl,nstl,DD,DDV,KERRE,INDIC,ITHHER)
 
      DO 13 I=1,nstl
        r_z = 0.D0
        DO J=1,nstl
          r_z = r_z + DD(I,J)*DEPST(J)
        ENDDO
        DSIGT(I)=r_z
 13   CONTINUE
 
C===============================================
C  CALCUL DE L INCREMENT DE DEFORMATIONS TOTALES
C===============================================
C
      CALL ELAST1(2,IFOUR,VAR0,NNVARI,XMAT1,NXMATT,YUNGV,XNUV,
     1          XCARB,ICARA,mfrl,nstl,DD,DDV,KERRE,INDIC,ITHHER)
      CALL ELAST1(2,IFOUR,VAR0,NNVARI,XMAT2,NXMATT,YUNGV,XNUV,
     1          XCARB,ICARA,mfrl,nstl,DDINV,DDV,KERRE,INDIC,ITHHER)
 
           AA=CTEPS
           DO 45 I=1,nstl
             IF (I.GT.3) AA=0.D0
             r_z = DEPST(I)
             DO 46 J=1,nstl
c*               r_z = r_z + (DDINV(I,J)*SIG0(J))-(DD(I,J)*SIG0(J))
               r_z = r_z + (DDINV(I,J)-DD(I,J))*SIG0(J)
46           CONTINUE
             DEPS(I)=r_z-AA
45         CONTINUE
C
C ------
C
      CALL INITTE(NYOG,YOG,NYNU,YNU,NYRHO,YRHO,NYALFA,YALFA,NYN,
     1     YN,NYM,YM,NYKK,YKK,NYKX,YKX,NNKX,NKX,X,T,
     2     XMAT,NXMATT,TRUC,NCOUR7)
 
 5    CONTINUE
      DO 10 I=1,nstl
                      SIGF(I)=SIG0(I)+DSIGT(I)
                      SIGMA(I)=SIG0(I)
   10 CONTINUE
      DO 20 I=1,NVARI
                      VARF(I)=VAR0(I)
   20 CONTINUE
C
      IF (ITHHER.EQ.2.AND.IALPH.GT.0) GOTO 32
   30 CONTINUE
       CALL CRIDAM(SIGF,nstl,VARF,XMAT,NXMATT,XCARB,ICARA,INPLAS,
     1NYKX,YKX,NNKX,NKX,NNVARI,mfrl,XF,KERRE,SIGGD,YY,T,TRUC,
     2                                                    NNFUS)
      IF (XF.LT.PRECI) GOTO 1000
      CALL CRIDOM(SIG0,nstl,VAR0,XMAT,NXMATT,XCARB,ICARA,INPLAS,
     1NNVARI,mfrl,XG,KERRE,SIGGD,Y)
C
C............ Y DESIGNE LE SEUIL D'ENDOMMAGEMENT .................
      IF (Y.EQ.ZEERO) THEN
                           INDIC=1
                           GOTO 32
      ENDIF
      TTEST=XG/Y
      IF (TTEST.GT.PRECIG) INDIC=1
 
C ==============
C ON A PLASTIFIE
C ==============
C ------------------------------------------------------------------
C DEFINITION DE DD, BB, R, S, XMU, ET QSI:
C  .          .         .
C SIG = DD*EPSELAS + BB*Q   AVEC:
C         - SIG     = CONTR.
C         - EPSELAS = DEFORM. ELASTIQUES
C         - DD      = MATR. DE HOOK ENDOMMAGEE
C         - Q       = VARIABLES INTERNES PILOTANT LES LOIS DU MODELE
C R   = LOI D'ECOULEMENT PLASTIQUE
C S   = LOIS D'EVOLUTION DES VARIABLES Q
C XMU = DERIVEE DU CRITERE PAR RAPPORT AUX CONTR.
C QSI = DERIVEE DU CRITERE PAR RAPPORT AUX VARIABLES Q
C ------------------------------------------------------------------
   32 ICO=0
C
C ------------------------------------------
C DEBUT DE LA BOUCLE DES ITERATIONS INTERNES
C --------------------------------------------------------------------
C A CHAQUE ITERATION, CORRESPONDANT A SIGF ET VARF CALCULES AU RANG I,
C "ENDOM1" CALCULE DD, BB, R, S, XMU ET QSI EN FONCTION DU MODELE DU
C MATERIAU ENDOMMAGEABLE.
C --------------------------------------------------------------------
   35 CONTINUE
 
      CALL TRACTT(YKX,NKX,TRUC,NYKX,NNKX,NNFUS,T)
      NFUST=2*NNFUS
      CALL DERTRA(NYOG,YOG,T,TYUNG,YUNGV,TINF,TSUP)
      CALL DERTRA(NYNU,YNU,T,TNU,XNUV,TINF,TSUP)
      CALL ENDOM1(SIGF,nstl,VARF,XMAT,NXMATT,XCARB,ICARA,
     1            NNFUS,TRUC,NNVARI,mfrl,IFOUR7,DD,DDINV,
     2            BB,R,S,XMU,QSI,KERRE,SIGGD,W,F,INDIC,ITHHER,
     3            YUNGV,XNUV)
      IF (KERRE.NE.0) GOTO 1000
C
      DO 50 J=1,nstl
        r_z = 0.D0
        DO 40 K=1,nstl
          r_z = r_z + DD(J,K)*R(K)
   40   CONTINUE
        DDR(J)= r_z
   50 CONTINUE
 
      DO 70 J=1,nstl
        r_z = 0.D0
        DO 60 K=1,NNVARI
          r_z = r_z + BB(J,K)*S(K)
   60   CONTINUE
        BBS(J) = r_z
   70 CONTINUE
 
      DLAMA1=0.D0
      DO 80 J=1,nstl
        DLAMA1=DLAMA1+XMU(J)*(DDR(J)-BBS(J))
   80 CONTINUE
      DLAMA2=0.D0
      DO 90 J=1,NNVARI
        DLAMA2=DLAMA2+QSI(J)*S(J)
   90 CONTINUE
 
C -------------------------------------------------------------
C CALCUL DE LA VARIATION DU MULT. PLAST. PAR LE CRIT. DE PLAST.
C -------------------------------------------------------------
C*?      r_DT=T-T0
      DLAMDA=XF/(DLAMA1-DLAMA2)
 
C ------------------------------------------------------------------
C CALCUL DES CONTR., DES VAR. INT. Q ET DE LA DEFORM. PLAST. CUMULEE
C AU RANG I+1
C ------------------------------------------------------------------
      DO 100 J=1,nstl
         SIGF(J)=SIGF(J)-DLAMDA*(DDR(J)-BBS(J))
  100 CONTINUE
      VARF(2)=VARF(2)+DLAMDA*S(1)
 
      IF (INDIC.EQ.1) THEN
      S2=S(2)
      IF (IFOUR7.NE.-2) THEN
      XA=1.D0+VARF(3)
      CALL ENDOM1(SIGF,nstl,VARF,XMAT,NXMATT,XCARB,ICARA,
     1            NNFUS,TRUC,NNVARI,mfrl,IFOUR7,DD,DDINV,
     2            BB,R,S,XMU,QSI,KERRE,SIGGD,W,F,INDIC,ITHHER,
     3            YUNGV,XNUV)
 
      IF (KERRE.NE.0) GOTO 1000
      XB=VARF(3)+DLAMDA*S(2)*(1.D0-VARF(3))
      XDELTA=(XA*XA)-4.D0*XB
            IF (XDELTA.LT.0.D0) THEN
               VARF(3)=VARF(3)+DLAMDA*S2
            ELSE
                XDELTA=SQRT(XDELTA)
                VARF(3)=(XA-XDELTA)/2.D0
            ENDIF
      ELSE
C
C  CAS DES CONTR. PLANES - MAJ DE D
C
      XA=VARF(3)+0.5D0*DLAMDA*S(2)
      CALL ENDOM1(SIGF,nstl,VARF,XMAT,NXMATT,XCARB,ICARA,
     1NNFUS,TRUC,NNVARI,mfrl,IFOUR7,DD,DDINV,
     2BB,R,S,XMU,QSI,KERRE,SIGGD,W,F,INDIC,ITHHER,
     3YUNGV,XNUV)
      XB=0.5D0*DLAMDA*S(2)*(1.D0-VARF(3))
      XDELTA=(1.D0-XA)*(1.D0-XA)-4.D0*XB
            IF (XDELTA.LT.0.D0) THEN
               VARF(3)=VARF(3)+DLAMDA*S2
            ELSE
               XDELTA=SQRT(XDELTA)
               VARF(3)=(1.D0+XA-XDELTA)/2.D0
            ENDIF
      ENDIF
      ENDIF
C
      VARF(1)=VARF(1)+DLAMDA/(1.D0-VARF(3))
C
C CONTR. PLANES - MAJ DES DEFORM. PLAST.
C
      IF (IFOUR7.EQ.-2) THEN
        CALL ENDOM1(SIGF,nstl,VARF,XMAT,NXMATT,XCARB,ICARA,
     1  NNFUS,TRUC,NNVARI,mfrl,IFOUR7,DD,DDINV,
     2  BB,R,S,XMU,QSI,KERRE,SIGGD,W,F,INDIC,ITHHER,
     3  YUNGV,XNUV)
        DO 110 J=1,nstl
              VARF(1+NNVARI+J)=VARF(1+NNVARI+J)+DLAMDA*R(J)
  110 CONTINUE
      ENDIF
 
C -------------------------------------------------------------
C SI LE CRIT. PLAST. EST ATTEINT, ON ARRETE LES ITERATIONS
C SI LE NBR. MAXIMUM D'ITERATIONS EST DEPASSE, IDEM MAIS ERREUR
C -------------------------------------------------------------
      CALL CRIDAM(SIGF,nstl,VARF,XMAT,NXMATT,XCARB,ICARA,INPLAS,
     1     NYKX,YKX,NNKX,NKX,NNVARI,mfrl,XF,KERRE,SIGGD,YY,T,TRUC,
     2                                                   NNFUS)
      TEST=ABS(XF/YY)
      IF (KERRE.NE.0) GOTO 1000
      IF (TEST.LT.PRECI.AND.INDIC.EQ.1) GOTO 120
      IF (TEST.LT.PRECI.AND.INDIC.EQ.0) THEN
      CALL CRIDOM(SIGF,nstl,VARF,XMAT,NXMATT,XCARB,ICARA,INPLAS,
     1            NNVARI,mfrl,XG,KERRE,SIGGD,Y)
             TTEST=XG/Y
             ATTEST=ABS(TTEST)
             IF (ATTEST.LT.PRECG) GOTO 120
             IF (TTEST.LE.PRECIG) THEN
                IF (ALPH.EQ.1.D0) GOTO 120
                A=1.D0
             ELSE
                A=-1.D0
             ENDIF
             ALPH=ALPH1+0.5D0*A*ABS(ALPH1-ALPH0)
             ALPH0=ALPH1
             ALPH1=ALPH
             IALPH=IALPH+1
             T=(1-ALPH)*T0+ALPH*T1
             X=VARF(2)
             CALL INITTE(NYOG,YOG,NYNU,YNU,NYRHO,YRHO,NYALFA,YALFA,
     1       NYN,YN,NYM,YM,NYKK,YKK,NYKX,YKX,NNKX,NKX,X,T,
     2       XMAT,NXMATT,TRUC,NCOUR7)
             ALFFA=XMAT(4)
C
             CALL ELAST1(1,IFOUR7,VAR0,NNVARI,XMAT,NXMATT,YUNGV,
     1       XNUV,XCARB,ICARA,mfrl,nstl,DDINV,DDV,KERRE,INDIC,ITHHER)
             CALL ELAST1(2,IFOUR7,VAR0,NNVARI,XMAT1,NXMATT,YUNGV,
     1       XNUV,XCARB,ICARA,mfrl,nstl,DD,DDV,KERRE,INDIC,ITHHER)
             CALL MULMAT(DDV,DDINV,DD,nstl,nstl,nstl)
C
             CTEPS=ALFA0*(T0-TREF)-ALFFA*(T-TREF)
             CALL ZDANUL(DSIGT,nstl)
C
             DO 8 I=1,nstl
             AA=1.D0
               DO 9 J=1,nstl
                IF (J.GT.3) AA=0.D0
                DSIGT(I)=DSIGT(I)+DDINV(I,J)*((ALPH*DEPS(J))+(AA*CTEPS))
                DSIGT(I)=DSIGT(I)+DDV(I,J)*SIG0(J)
9              CONTINUE
               DSIGT(I)=DSIGT(I)-SIG0(I)
8            CONTINUE
             IF (IALPH.GT.50) GOTO 1000
             GOTO 5
      ENDIF
      ICO=ICO+1
      IF (ICO.GT.ITMAX) THEN
                            KERRE=2
                            GOTO 1000
      ENDIF
      GOTO 35
C ----------------------------------------
C FIN DE LA BOUCLE DES ITERATIONS INTERNES
C ----------------------------------------
C
C ====================================================
C FIN DES ITERATIONS INTERNES - LE CRITERE EST VERIFIE
C                 MAJ DES VARIABLES
C ====================================================
C
C ----------------------------------------------------------------
C CALCUL DES DEFORM. ELAST. A LA FIN DU PAS - ELLES SERONT DS. DDR
C ----------------------------------------------------------------
  120 CONTINUE
      CALL ELAST1(2,IFOUR7,VARF,NNVARI,XMAT,NXMATT,YUNGV,XNUV,
     &       XCARB,ICARA,mfrl,nstl,DD,DDV,KERRE,INDIC,ITHHER)
      DO 140 J=1,nstl
        r_z = 0.D0
        DO 130 K=1,nstl
          r_z = r_z + DD(J,K)*SIGF(K)
  130   CONTINUE
        DDR(J) = r_z
  140 CONTINUE
C
C -----------------------------------------
C CALCUL DES DEFORM. PLAST. A LA FIN DU PAS
C ----------------------------------------------------------------
C ON CALCULE D'ABORD, DS. BBS, LES DEFORM. TOTALES A LA FIN DU PAS
C ----------------------------------------------------------------
      CALL ELAST1(2,IFOUR7,VAR0,NNVARI,XMAT1,NXMATT,YUNGV,XNUV,
     1XCARB,ICARA,mfrl,nstl,DD,DDV,KERRE,INDIC,ITHHER)
      CALL ZDANUL(BBS,nstl)
      DO 170 J=1,nstl
      DO 160 K=1,nstl
                      BBS(J)=BBS(J)+DD(J,K)*SIGMA(K)
  160 CONTINUE
  170 CONTINUE
      AA=CTEPS
      DO 180 J=1,nstl
      IF (J.GT.3) AA=0.D0
      BBS(J)=BBS(J)+VAR0(1+NNVARI+J)+(ALPH*DEPS(J))+(AA)
  180 CONTINUE
C
C
      DO 190 J=1,nstl
               IF (IFOUR7.EQ.-2) GOTO 195
               VARF(1+NNVARI+J)=BBS(J)-DDR(J)
  195          DEFP(J)=VARF(1+NNVARI+J)-VAR0(1+NNVARI+J)
  190 CONTINUE
      IF (ALPH.EQ.1.D0) GOTO 900
      IF (INDIC.EQ.0) THEN
                          INDIC=1
                          ALPH=1.D0-ALPH
                         DO 520 J=1,NVARI
                          VAR0(J)=VARF(J)
  520                    CONTINUE
                         DO 525 J=1,nstl
                          SIGMA(J)=SIGF(J)
  525                    CONTINUE
C
                          T0=T
                          T=T1
                          DO 530 J=1,NXMATT
                          XMAT1(J)=XMAT(J)
                          XMAT(J)=XMAT2(J)
  530                     CONTINUE
                          ALFA0=XMAT1(4)
                          ALFFA=XMAT(4)
             CALL ELAST1(1,IFOUR7,VAR0,NNVARI,XMAT,NXMATT,YUNGV,
     1       XNUV,XCARB,ICARA,mfrl,nstl,DD,DDV,KERRE,INDIC,ITHHER)
             CALL ELAST1(2,IFOUR7,VAR0,NNVARI,XMAT1,NXMATT,YUNGV,
     1       XNUV,XCARB,ICARA,mfrl,nstl,DDINV,DDV,KERRE,INDIC,
     2       ITHHER)
             CALL MULMAT(DDV,DD,DDINV,nstl,nstl,nstl)
C
             CTEPS=ALFA0*(T0-TREF)-ALFFA*(T-TREF)
C
             DO 515 I=1,nstl
               AA=CTEPS
               r_z = 0.D0
               DO 510 J=1,nstl
                 IF (J.GT.3) AA=0.D0
                 r_z = r_z + DD(I,J)*((ALPH*DEPS(J))+AA)
     &                     + DDV(I,J)*SIGMA(J)
510            CONTINUE
               DSIGT(I)=r_z-SIGMA(I)
               SIGF(I)=SIGF(I)+DSIGT(I)
515          CONTINUE
                          GOTO 30
      ENDIF
C
C ==================================================================
C SI LE PT. DE GAUSS EST ROMPU, ON ANNULE PROGRESSIVEMENT SES CONTR.
C ==================================================================
  900 DCT=XMAT(7)
      DCC=1.1D0*DCT
      XD=VARF(3)
      IF (XD.GE.DCC) THEN
         VARF(3)=1.D0
         CALL ZDANUL(SIGF,nstl)
      ENDIF
      IF (XD.GT.DCT.AND.XD.LT.DCC) THEN
         VARF(3)=1.D0
         BETA=11.D0-10.D0*(XD/DCT)
         DO J=1,nstl
            SIGF(J)=SIGF(J)*BETA
         ENDDO
 
      ENDIF
 1000 RETURN
      END
*******************************************************************