Sources Esope | Cast3M

cconst
C CCONST    SOURCE    PV        22/04/22    21:15:05     11344          
      SUBROUTINE CCONST(wrk52,wrk53,wrk54,WRK7,WRK8,WRK9,
     1     WRK91,NVARI,NSSINC,INV,IFOURB,T0,TF,FI0,FIF
     4    ,TLIFE,NCOURB,IB,IGAU,NBPGAU,KERREU1,iecou,xecou)
C CONSTI    SOURCE    BROC      00/12/20    21:15:56     4058
c      SUBROUTINE CONSTI(WRK0,WR00,WRK1,WRK5,WRK7,WRK8,WRK9,WTRAV,
c     1 INPLAS,MFR,DT,NSTRS,NVARI,NCOMAT,PRECAS,MSOUPA,JECHER,DTT,
c     2 NSSINC,INV,KERRE,ICARA,IFOURB,NYOG,NYNU,NYALFA,NYSMAX,
c     3 NYN,NYM,NYKK,NYALF1,NYBET1,NYR,NYA,NYKX,NNKX,NYRHO,NSIGY,T0,TF,
c     5 TREF,TLIFE,ITHHER,NCOURB,CMATE,N2EL,N2PTEL,IB,IGAU,EPAIST,
c     7NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,KERREU1)
C
C---------------------------------------------------------------------
C Objet: Calculer au cours d'un pas de temps DT, l'evolution des
C        variables internes a l'aide d'un schema Runge-Kutta 1.2
C ---------------------------------------------------------------------
C MFR1 &lt;- MFR, XCARB &lt;- XCAR, NSTRS1 &lt;- NSTRS,
C
C---------------------------------------------------------------------
C Entree: INPLAS type de materiau
C         MFR indice de la formulation mecanique(seulement massif ou coque
C           pour les materiaux endommageables)
C         DEPST(NSTRS1) increment des deformations totales
C         SIG0(NSTRS1) contraintes initiales
C         EPIN0(NSTRS1) deformations viscoplastiques initiales
C         VAR0(NVARI) variables internes initiales
C         NVARI nombre de variables internes
C         YOG(NYOG) courbe du module d'Young en fonction de T°C
C         YNU(NYNU) courbe du coefficient de Poisson en fonction de T°C
C         SIGY(NSIGY) courbe de la limite elastique en fonction de T°C
C         YRHO(NYRHO) courbe de la masse volumique en fonction de T°C
C         YALFA(NYALFA) courbe du coeff de dilatation en fonction de T°C
C         YN(NYN)
C         YM(NYM)
C         YKK(NYKK)
C         YALFA1(NYALF1)   courbes des autres coefficients caracteristiques
C         YBETA1(NYBET1)   en fonction de la T°C intervenant
C         YALF2(NYALF2)    dans les lois d'evolution
C         YBET2(NYBET2)
C         YR(NYR)
C         YA(NYA)
C         YKX(NYKX) fonction k(X) tabulee en fonction de la temperature
C         NKX(NNKX) tableau de pointeurs sur les courbes de k(X)
C         XMAT(NCOMAT) materiau
C         XCARB(ICARA) caracteristiques geometriques
C         YSMAX(NYSMAX) intervient ds. le test de convergence des iter.
C         TRUC(NCOURB) tableau de travail
C         PRECIS precision relative sur SIGMA
C         MSOUPA nombre maximal de sous pas autorises
C         JECHER = 0 avancer
C                = 1 rechercher sortie avec DTT
C         IFOURB = -3 EN DEFORM. PLANES GENER.
C                  -2 EN CONTR.PLANES
C                  -1 EN DEFORM. PLANES
C                   0 EN AXISYMETRIE
C                   1 EN SERIE DE FOURIER
C                   2 EN TRIDIM
*         CMATE = NOM DU MATERIAU
*         VALMAT= TABLEAU DE CARACTERISTIQUES DU MATERIAU
*         VALCAR= TABLEAU DE CARACTERISTIQUES GEOMETRIQUES
*         N2EL  = NBRE D ELEMENTS DANS SEGMENT DE HOOKE
*         N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE
*         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         DD(NSTRS1,NSTRS1) matrice de Hooke en fonction de T
C         DDV(NSTRS1,NSTRS1) derivee de DD
C         DDINV(NSTRS1,NSTRS1) inverse de DD
C         T0 temperature a t0
C         TF temperature a t0+DT
C         FI0 densite de fissions a t0
C         FIF densite de fissions a t0+DT
C         TREF temperature de reference
C         DT pas de temps
C         ITHHER = 0 pas de chargement thermique et materiau constant
C                = 1 chargement thermique et materiau constant
C                = 2 chargement thermique et materiau(T)
C-----------------------------------------------------------------------
C
C-----------------------------------------------------------------------
C Sortie: SIGF(NSTRS1) contraintes finales
C         EPINF(NSTRS1) deformations viscoplastiques finales
C         VARF(NVARI) variables internes finales
C         DTT sous-increment de temps optimal (si JECHER=1)
C         TLIFE  sous-increment de temps a rupture pour materiau
C           viscoplastique endommageable
C         NSSINC nombre de sous-increments si JECHER=0
C         INV = 1 si inversion
C               0 sinon
C         KERRE = 0 si tout OK
C               &lt;> 0 si entrees incoherentes
C-----------------------------------------------------------------------
C
      IMPLICIT INTEGER(I-N)
       IMPLICIT REAL*8(A-H,O-Z)
-INC CCREEL
 
-INC PPARAM
-INC CCOPTIO
-INC DECHE
 
 
      SEGMENT IECOU
*      COMMON/IECOU/NYOG,NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK,
       INTEGER     NYOG, NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK,
C      INTEGER NYOG, NYNU, NYALFA,NYSMAX,NYN, NYM,  NYKK,
     1  NYALF1,NYBET1,NYR , NYA,   NYRHO,NSIGY, NNKX,NYKX,icow16,
C    .  NYALF1,NYBET1,NYR,  NYA,   NYRHO,NSIGY, NNKX,NYKX,  IND,
     2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24,
C    .  NSOM,  NINV, NINCMA,NCOMP, JELEM, LEGAUS,INAT,  NCXMAT,
     3 icow25,icow26,icow27,icow28,icow29,icow30,ICARA,
C    .  LTRAC, MFR,  IELE,  NHRM,   NBNN, NBELEM,ICARA,
     4  icow32,icow33,NSTRS1,MFR1,  NBGMAT,NELMAT,MSOUPA,
C    .  LW2,   NDEF,  NSTRSS,MFR1,  NBGMAT,NELMAT,MSOUPA,
     5  icow39,icow40,icow41,icow42,icow43,icow44,NYOG1,
C    .  NUMAT1,LENDO, NBBB,  NNVARI,KERR1, MELEME,NYOG1,
     6  NYNU1,NYALFT1,NYSMAX1,NYN1,NYM1,NYKK1,NYALF2,
C    .  NYNU1,NYALFT1,NYSMAX1,NYN1,NYM1,NYKK1,NYALF2,
     7  NYBET2,NYR1,NYA1,NYQ1,NYRHO1,NSIGY1
C    .  NYBET2,NYR1,NYA1,NYQ1,NYRHO1,NSIGY1
      ENDSEGMENT
 
 
      SEGMENT XECOU
*      COMMON/XECOU/DTOPTI,TSOM,TCAR,DTT,DT,TREFA,TEMP00
      REAL*8 xcow1, xcow2,xcow3,DTT  ,DT, TREF, xcow7
C     REAL*8 DTOPTI,TSOM, TCAR, DTT,  DT,   TREFA,TEMP00
      ENDSEGMENT
 
 
C
C
      SEGMENT WRK7
        REAL*8 F(NCOURB,2),W(NCOURB),TRUC(NCOURB)
      ENDSEGMENT
C
      SEGMENT WRK8
        REAL*8 DD(NSTRS1,NSTRS1),DDV(NSTRS1,NSTRS1),DDINV(NSTRS1,NSTRS1)
        REAL*8 DDINVp(NSTRS1,NSTRS1)
      ENDSEGMENT
C
      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
C
      SEGMENT WRK91
        REAL*8 YOG1(NYOG1),YNU1(NYNU1),YALFT1(NYALFT1),YSMAX1(NYSMAX1)
        REAL*8 YN1(NYN1),YM1(NYM1),YKK1(NYKK1),YALF2(NYALF2)
        REAL*8 YBET2(NYBET2),YR1(NYR1),YA1(NYA1),YQ1(NYQ1),YRHO1(NYRHO1)
        REAL*8 SIGY1(NSIGY1)
      ENDSEGMENT
C
       DIMENSION VAR(100),VAR1(100),VARP1(100),VARP2(100)
       DIMENSION CRIGI(12),VAR12(100),VARP3(100)
       DIMENSION VARP4(100),VAR13(100),VART(100)
       DIMENSION VART1(100),ZMAT(100)
C
C#MC 21/01/99 : les tableaux doivent dimensionnes en fonction
C               du plus grand INFELE(16) (voir elquoi.eso)
       DIMENSION SIG(8),SIG1(8),SIG12(8),SIG13(8)
       DIMENSION DSPT(8),EPSTHD(8),XX(8)
       DIMENSION EVP1(8),EVP2(8),XPM1(8),XPM2(8),EVP3(8),EVP4(8)
       DIMENSION XPM3(8),XPM4(8)
       DIMENSION SIGP1(8),SIGP2(8), SIGP3(8), SIGP4(8)
       DIMENSION EPSV(8),EPSV1(8),EPSV12(8),EPSV13(8)
C
       logical dtlibr,iforce,ilog
C
       CALL ZERO(SIG1 ,8    ,1)
C
       NCOMAT = NMATT
C
       TPS0 = TEMP0
       TPSF = TEMPF
       TPSX = MAX(ABS(TPS0),ABS(TPSF))
       DT0  = TPSF - TPS0
       IF (ABS(DT0).LE.ABS(XZPREC*TPSX)) DT0 = 0.D0
 
       iffo=0
       IVTEST=0
*       pasbea=0.d0
       PRELOC=1.d-8
       msoupa=1000000
*      write(6,*) 'cconst ncomat',ncomat
*      write(6,*) (xmat(iu),iu=1,ncomat+1)
        dtlibr=.TRUE.
C Test sur l'identite de toutes les listes de temperatures des coefficients
C intervenant dans les lois d'evolutions non-lineaires des variables internes
       IF (INPLAS.EQ.29) THEN
*       write(6,*) 'cconst avant test kerre',kerre
         CALL TEST(YN,NYN,YM,NYM,YKK,NYKK,YALFA1,NYALF1,YBETA1,NYBET1,
     &             YR,NYR,YA,NYA,SIGY,NSIGY,YKX,NYKX,NKX,NNKX,ITEST)
*     write(6,*) 'cconst avant test kerre',kerre
       ENDIF
       IF (INPLAS.EQ.142) THEN
*        write(6,*) 'cconst avant test1 kerre',kerre
         CALL TEST1(YN1,NYN1,YM1,NYM1,YKK1,NYKK1,YALF2,NYALF2,YBET2,
     &   NYBET2,YR1,NYR1,YA1,NYA1,SIGY1,NSIGY1,YQ1,NYQ1,ITEST)
*        write(6,*) 'cconst avant test1 kerre',kerre
       ENDIF
C
       KERRE = 0
       IF (MFR1.NE.1.AND.MFR1.NE.3.AND.MFR1.NE.5.AND.MFR1.NE.17.AND.
     &     MFR1.NE.31.AND.MFR.NE.33) THEN
         KERRE = 99
         RETURN
       ENDIF
*
*  AM 5/5/00   MFR1 = 33   : MODELES 19 A 24 pour le moment
*  AM 27/5/3               :  ET  44,45
*
       IF (MFR1.EQ.33) THEN
         IF ((INPLAS.LT.19.OR.INPLAS.GT.24).AND.
     &       (INPLAS.NE.44.AND.INPLAS.NE.45))THEN
           KERRE = 99
           RETURN
         ENDIF
       ENDIF
*
*      write(6,*) 'cconst mfr1',mfr1
       IF (MFR1.EQ.3) THEN
         THICK = XCARB(1)
         ALFA  = XCARB(2)
       ENDIF
*       dtprem=0.D0
*       dtdeux=0.d0
       dtleft= dt
       BORNE = 2.0
       RMAX  = 1.3
       RMIN  = 0.7
       DIV   = 7.0
       FAC   = 3.0
       TLIFE = -1.D0
C
C   CALCUL DES INCREMENTS DE DEFORMATIONS
C
       IF ((INPLAS.NE.29).AND.(INPLAS.NE.142)) THEN
         CALL CALSIG(DEPST,DDAUX,NSTRS1,CMATE,VALMAT,VALCAR,N2EL,
     1               N2PTEL,MFR1,IFOURB,IB,IGAU,EPAIST,NBPGAU,MELE,
     2               NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,XLOC,XGLOB,
     3               D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD)
C
         IF (IRTD.NE.1) THEN
           KERRE=69
           GOTO 998
         ENDIF
C
       ENDIF
       IF (MFR1.EQ.3) THEN
         DO 10 I=1,NSTRS1/2
           SIG0(         I)= SIG0(         I)/THICK
           SIG0(NSTRS1/2+I)= SIG0(NSTRS1/2+I)*6.0D0/THICK/THICK
           IF ((INPLAS.NE.29).AND.(INPLAS.NE.142)) THEN
             DSIGT(         I)= DSIGT(         I)/THICK
             DSIGT(NSTRS1/2+I)= DSIGT(NSTRS1/2+I)*6.0D0/THICK/THICK
           ELSE
             DEPST(         I)=  DEPST(         I)
             DEPST(NSTRS1/2+I)= -DEPST(NSTRS1/2+I)*THICK/2.D0
           ENDIF
10       CONTINUE
         IF (IFOURB.EQ.-2) THEN
           SIG0(2)=0.D0
           SIG0(4)=0.D0
           IF ((INPLAS.NE.29).AND.(INPLAS.NE.142)) THEN
             DSIGT(2)=0.D0
             DSIGT(4)=0.D0
           ENDIF
         ENDIF
       ENDIF
C
C REMISE A ZERO DE SIG A L'EXCEPTION DU MOMENT SUIVANT Z ET DE
C   L'EFFORT SUIVANT X (MODE I DU CHARGEMENT)
C
       IF (MFR1.EQ.17) THEN
         SIG0(2) = 0.D0
         SIG0(3) = 0.D0
         SIG0(4) = 0.D0
         SIG0(5) = 0.D0
       ENDIF
C
C------------------------------------------
C CONTROLE DE LA COHERENCE DES ENTREES
C------------------------------------------
       IF (DT.LT.0.0) KERRE = 414
       IF (INPLAS.EQ.63.AND.MFR1.NE.1.AND.MFR1.NE.31) THEN
         KERRE=99
         RETURN
       ENDIF
       IF (DT.EQ.0.0) DT = 1.e-20
       MOTERR(1:8) = ' CONST  '
       IF (INPLAS.EQ.17) THEN
         IF ((NVARI.NE.(6+4*NSTRS1)).AND.(MFR1.NE.5)) KERRE = 146
         IF ((NVARI.NE.(10+4*NSTRS1)).AND.(MFR1.EQ.5)) KERRE = 146
         IF (IFOURB.NE.-2.AND.NCOMAT.LT.24) KERRE = 146
         IF (IFOURB.EQ.-2.AND.NCOMAT.LT.25) KERRE = 146
         XMAX=XMAT(8)
         GOTO 30
       ENDIF
 
       IF (MFR1.NE.33) THEN
*           write(6,*) 'cconst inplas ifourb ncomat',inplas,ifourb,ncomat
         IF (INPLAS.EQ.19.AND.IFOURB.NE.-2.AND.NCOMAT.LT. 8)KERRE = 146
         IF (INPLAS.EQ.19.AND.IFOURB.EQ.-2.AND.NCOMAT.LT. 9)KERRE = 146
         IF (INPLAS.EQ.20.AND.IFOURB.NE.-2.AND.NCOMAT.LT.18)KERRE = 146
         IF (INPLAS.EQ.20.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.19)KERRE = 146
         IF (INPLAS.EQ.21.AND.IFOURB.NE.-2.AND.NCOMAT.LT.12)KERRE = 146
         IF (INPLAS.EQ.21.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.13)KERRE = 146
         IF (INPLAS.EQ.22.AND.IFOURB.NE.-2.AND.NCOMAT.LT.12)KERRE = 146
         IF (INPLAS.EQ.22.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.14)KERRE = 146
         IF (INPLAS.EQ.23.AND.IFOURB.NE.-2.AND.NCOMAT.LT.16)KERRE = 146
         IF (INPLAS.EQ.23.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.17)KERRE = 146
         IF (INPLAS.EQ.24.AND.IFOURB.NE.-2.AND.NCOMAT.LT.11)KERRE = 146
         IF (INPLAS.EQ.24.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.12)KERRE = 146
         IF (INPLAS.EQ.25.AND.IFOURB.NE.-2.AND.NCOMAT.LT.26)KERRE = 146
         IF (INPLAS.EQ.25.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.27)KERRE = 146
C         IF (INPLAS.EQ.29.AND.IFOURB.NE.-2.AND.NCOMAT.LT.13)KERRE = 146
C         IF (INPLAS.EQ.29.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.14)KERRE = 146
         IF (INPLAS.EQ.44.AND.IFOURB.NE.-2.AND.NCOMAT.LT.20)KERRE = 146
         IF (INPLAS.EQ.44.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.21)KERRE = 146
         IF (INPLAS.EQ.45.AND.IFOURB.NE.-2.AND.NCOMAT.LT.27)KERRE = 146
         IF (INPLAS.EQ.45.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.28)KERRE = 146
         IF (INPLAS.EQ.53.AND.IFOURB.NE.-2.AND.NCOMAT.LT.28)KERRE = 146
         IF (INPLAS.EQ.53.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.29)KERRE = 146
         IF (INPLAS.EQ.61.AND.IFOURB.NE.-2.AND.NCOMAT.LT.18)KERRE = 146
         IF (INPLAS.EQ.61.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.19)KERRE = 146
         IF (INPLAS.EQ.63.AND.IFOURB.NE.-2.AND.NCOMAT.LT.32)KERRE = 146
         IF (INPLAS.EQ.63.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.33)KERRE = 146
         IF (INPLAS.EQ.70.AND.IFOURB.NE.-2.AND.NCOMAT.LT.14)KERRE = 146
         IF (INPLAS.EQ.70.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.15)KERRE = 146
         IF (INPLAS.EQ.76.AND.IFOURB.NE.-2.AND.NCOMAT.LT.26)KERRE = 146
         IF (INPLAS.EQ.76.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.27)KERRE = 146
         IF (INPLAS.EQ.77.AND.IFOURB.NE.-2.AND.NCOMAT.LT.18)KERRE = 146
         IF (INPLAS.EQ.77.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.19)KERRE = 146
         IF (INPLAS.EQ.83.AND.IFOURB.NE.-2.AND.NCOMAT.LT.15)KERRE = 146
         IF (INPLAS.EQ.83.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.17)KERRE = 146
         IF (INPLAS.EQ.84.AND.IFOURB.NE.-2.AND.NCOMAT.LT.13)KERRE = 146
         IF (INPLAS.EQ.84.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.14)KERRE = 146
         IF (INPLAS.EQ.85.AND.IFOURB.NE.-2.AND.NCOMAT.LT.19)KERRE = 146
         IF (INPLAS.EQ.85.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.20)KERRE = 146
         IF (INPLAS.EQ.86.AND.IFOURB.NE.-2.AND.NCOMAT.LT.17)KERRE = 146
         IF (INPLAS.EQ.86.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.18)KERRE = 146
         IF (INPLAS.EQ.102.AND.IFOURB.NE.-2.AND.NCOMAT.LT.25)KERRE = 146
         IF (INPLAS.EQ.102.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.26)KERRE = 146
         IF (INPLAS.EQ.130.AND.IFOURB.NE.-2.AND.NCOMAT.LT.10)KERRE = 146
         IF (INPLAS.EQ.130.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.11)KERRE = 146
         IF (INPLAS.EQ.136.AND.IFOURB.NE.-2.AND.NCOMAT.LT.10)KERRE = 146
         IF (INPLAS.EQ.136.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.11)KERRE = 146
         IF (INPLAS.EQ.137.AND.IFOURB.NE.-2.AND.NCOMAT.LT.10)KERRE = 146
         IF (INPLAS.EQ.137.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.11)KERRE = 146
         IF (INPLAS.EQ.138.AND.IFOURB.NE.-2.AND.NCOMAT.LT.10)KERRE = 146
         IF (INPLAS.EQ.138.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.11)KERRE = 146
         IF (INPLAS.EQ.139.AND.IFOURB.NE.-2.AND.NCOMAT.LT.10)KERRE = 146
         IF (INPLAS.EQ.139.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.11)KERRE = 146
 
C
       ELSE
C
C        cas MFR1=33
C
 
         IF (INPLAS.EQ.19.AND.IFOURB.NE.-2.AND.NCOMAT.LT.16)KERRE = 146
         IF (INPLAS.EQ.19.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.17)KERRE = 146
         IF (INPLAS.EQ.20.AND.IFOURB.NE.-2.AND.NCOMAT.LT.26)KERRE = 146
         IF (INPLAS.EQ.20.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.27)KERRE = 146
         IF (INPLAS.EQ.21.AND.IFOURB.NE.-2.AND.NCOMAT.LT.20)KERRE = 146
         IF (INPLAS.EQ.21.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.21)KERRE = 146
         IF (INPLAS.EQ.22.AND.IFOURB.NE.-2.AND.NCOMAT.LT.20)KERRE = 146
         IF (INPLAS.EQ.22.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.22)KERRE = 146
         IF (INPLAS.EQ.23.AND.IFOURB.NE.-2.AND.NCOMAT.LT.24)KERRE = 146
         IF (INPLAS.EQ.23.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.25)KERRE = 146
         IF (INPLAS.EQ.24.AND.IFOURB.NE.-2.AND.NCOMAT.LT.19)KERRE = 146
         IF (INPLAS.EQ.24.AND.IFOURB.EQ.-2.AND.NCOMAT.LT.20)KERRE = 146
C
       ENDIF
C
       IF (IFOURB.EQ.1) THEN
         KERRE = 194
         MOTERR(1:8) = 'FLUAGE'
       ENDIF
       XMAX=XMAT(5)
       IF ((INPLAS.EQ.25).OR.(INPLAS.EQ.53)) XMAX=XMAT(7)
       IF ((INPLAS.EQ.76).OR.(INPLAS.EQ.77)) XMAX=XMAT(7)
       IF ((INPLAS.EQ.70).OR.(INPLAS.EQ.107)) XMAX=XMAT(1)*1.D-3
       IF (INPLAS.EQ.29) THEN
         CALL DERTRA(NYSMAX,YSMAX,TF,XMAX,XMAXV,TO,TO)
       ENDIF
       IF (INPLAS.EQ.142) THEN
         CALL DERTRA(NYSMAX1,YSMAX1,TF,XMAX,XMAXV,TO,TO)
       ENDIF
C
C  TEST SUR XMAX          MILL  8/3/91
C
       IF (XMAX.EQ.0.D0) THEN
         IF (INPLAS.EQ.29) THEN
           CALL DERTRA(NYOG,YOG,TF,XMAX,XMAXV,TO,TO)
           XMAX=XMAX*1.D-3
         ELSEIF (INPLAS.EQ.142) THEN
           CALL DERTRA(NYOG1,YOG1,TF,XMAX,XMAXV,TO,TO)
           XMAX=XMAX*1.D-3
         ELSEIF ((INPLAS.NE.29).AND.(INPLAS.NE.142)) THEN
           XMAX=XMAT(1)*1.D-3
         ENDIF
       ENDIF
C
30     CONTINUE
C
C-----------------------------
       IF (KERRE.NE.0) THEN
         GOTO 999
       ENDIF
 
C
C===========================================================
C  A PARTIR DE MAINTENANT, LES DEFORMATIONS
C  DE CISAILLEMENT NE SONT PLUS
C  DEFINIES PAR DES GAMA.
C  ON DIVISE DONC LES TERMES DE CISAILLEMENT PAR 2.
C  CECI NE CONCERNE PAS LE MODELE VISCO-ENDOMMAGEABLE
C  DE LEMAITRE (INPLAS=29 ET INPLAS=142).
C
C  SEULES LES FORMULATIONS SUIVANTES SONT ACCEPTEES PAR CONSTI:
C      MFR1=1  (MASSIF)
C      MFR1=5  (COQUES EPAISSES)
C      MFR1=3  (COQUES MINCES)
C      MFR1=17 (TUYAUX FISSURES)
C      MFR1=31 (BBAR)
C      MFR1=33 (POREUX)
C
       IF ((INPLAS.NE.29).AND.(INPLAS.NE.142)) THEN
C
C  Cas de la formulation massive
C  Les termes de cisaillement apparaissent
C  au delà de la troisieme composante
C
         IF (MFR1.EQ.1.OR.MFR1.EQ.31.OR.MFR1.EQ.33) THEN
           DO I=4,NSTRS1
             EPIN0(I)=0.5D0*EPIN0(I)
           ENDDO
C
C  Cas des coques épaisses
C  Les termes de cisaillement apparaissent
C  au delà de la deuxieme composante
C
         ELSE IF (MFR1.EQ.5) THEN
           DO I=3,NSTRS1
             EPIN0(I)=0.5D0*EPIN0(I)
           ENDDO
C
C  Cas des coques minces
C  Les termes de cisaillement apparaissent
C  pour la troisieme et la sixieme composante
C  uniquement dans les cas de calculs
C  tridimensionnels ou d'analyse de Fourier
C
         ELSE IF (MFR1.EQ.3) THEN
           IF ((IFOURB.EQ.1).OR.(IFOURB.EQ.2)) THEN
             EPIN0(3)=0.5D0*EPIN0(3)
             EPIN0(6)=0.5D0*EPIN0(6)
           ENDIF
C
C  Reste le cas des tuyaux fissurés (MFR1=17)
C
         ENDIF
       ENDIF
C
C===========================================================
C
C ----------------
C  INITIALISATION
C ----------------
       ITERO = 0
 6543  CONTINUE
 
       itero = 1 + itero
       if ( itero.ne.1) THEN
           if(ib.eq.1.and.igau.eq.1) write(6,*)  'itero ', itero
         dtlibr = .true.
         preloc = preloc * 7.d0
c          write(6,*)  ' precision modifiée ', preloc
         if (itero.gt.3) then
****          kerre = 460
           kerre = 268
           return
         endif
       endif
       DTLEFT = DT
       TAU    = DTLEFT
       dtaumi= dtleft / 1500.
       TI0=T0
       TI1=TF
       TPOINT=(TF-T0)/DT
       FII0=FI0
       FII1=FIF
       FPOINT=(FIF-FI0)/DT
       ASIG = SQRT(PROCON(SIG0,SIG0,NSTRS1))
       ERRABS = PRELOC*ASIG
       IF (XMAX.GT.ASIG) ERRABS = PRELOC*XMAX
C SP   write(6,*) ' ### CCONST : Initialisation de SIG et EPSV'
C SP   write(6,*) ' ### CCONST : SIG0(1) =', SIG0(1)
       DO 40 I=1,NSTRS1
         SIG(I) = SIG0(I)
         EPSV(I) = EPIN0(I)
         IF ((INPLAS.NE.29).AND.(INPLAS.NE.142)) THEN
           DSPT(I) = DSIGT(I)/DT
         ENDIF
40     CONTINUE
*
*          iter=0
*       if ( .NOT.DTLIBR) THEN
*           dtminl = ( dt * 1.001 ) / msopua
*           r = dtseco / dtprem
* 1245      continue
*          iter=iter+1
**          if ( abs ( 1.- R) . gt. 0.001 ) then
*           bb = ( 1 - r**(msoupa-10)) / ( 1 - r) * dtprem
*          else
*           r = 1.d0
*           bb = dtprem * ( msoupa - 10)
*          endif
*          tau = dtprem * dt / bb * 1.0001
*          write (6,1234)iter,dtprem,dtseco,dtdeux,r,bb,dt,tau
*          if (bb . lt . dt/1.2) THEN
*              if ( iter.lt.15) then
*                r = r + abs ( 1. - r ) / 10.
*              else
*               kerre = 460
*               return
*             endif
*             go to 1245
*         endif
* 1234 format ('it pr se de r b t ta',i2, 7e9.3)
*
*           dtx = dt * 1.00001;
*           write(6,*) 'avpremdeux r', dtprem,dtdeux,r
*          call decoup(-msoupa+10,dtprem/dtx,dtdeux/dtx,r,nn,xde
*     $          ,xdf,dtx)
*           write(6,*) 'xde xdf de r',xde,xdf,r
*          tau = xde / r
*
*       endif
C
       IF (INPLAS.EQ.29) THEN
C
C================================================
C Calcul de l increment de deformation totale reel.
C On enleve donc tous les termes qui correspondent
C a l influence de la temperature et de
C l endommagement (travail inverse de ce qui est
C fait dans le procedure increme).
C================================================
C
********* materiau dependant de la temperature ***********************
         CALL DERTRA(NYOG,YOG,T0,YUNG0,YUNGV0,TO,TO)
         CALL DERTRA(NYNU,YNU,T0,ENU0,ENUV0,TO,TO)
         XMAT(1)=YUNG0
         XMAT(2)=ENU0
C------------------------------------------------
C Calcul de la matrice de Hooke inverse DD a t=t0
C------------------------------------------------
         CALL ELAST1(2,IFOURB,VAR0,NVARI,XMAT,NCOMAT,YUNGV0,ENUV0,
     &               XCARB,ICARA,MFR1,NSTRS1,DD,DDV,KERRE,2,ITHHER)
C
********* materiau dependant de la temperature ***********************
         CALL DERTRA(NYOG,YOG,TF,YUNG1,YUNGV1,TO,TO)
         CALL DERTRA(NYNU,YNU,TF,ENU1,ENUV1,TO,TO)
         XMAT(1)=YUNG1
         XMAT(2)=ENU1
C---------------------------------------------------
C Calcul de la matrice de Hooke inverse DDINV a T=TF
C---------------------------------------------------
         CALL ELAST1(2,IFOURB,VAR0,NVARI,XMAT,NCOMAT,YUNGV1,ENUV1,
     &               XCARB,ICARA,MFR1,NSTRS1,DDINV,DDV,KERRE,2,ITHHER)
***********************************************************************
C
         CALL DERTRA(NYALFA,YALFA,T0,ALFA0,ALFAV0,TO,TO)
         CALL DERTRA(NYALFA,YALFA,TF,ALFAF,ALFAVF,TO,TO)
         CTEPS=ALFA0*(T0-TREF)-ALFAF*(TF-TREF)
C
         CALL ZDANUL(DSPT,NSTRS1)
         AA=1.D0
         DO 45 I=1,NSTRS1
           DSPT(I)=DEPST(I)
           IF (I.GT.3) AA=0.D0
           DO 46 J=1,NSTRS1
             DSPT(I)=DSPT(I)+(DDINV(I,J)*SIG0(J))
             DSPT(I)=DSPT(I)-(DD(I,J)*SIG0(J))
46         CONTINUE
           DSPT(I)=DSPT(I)-(AA*CTEPS)
           DSPT(I)=DSPT(I)/DT
45       CONTINUE
C
       ELSEIF (INPLAS.EQ.142) THEN
C
C================================================
C Calcul de l increment de deformation totale reel.
C On enleve donc tous les termes qui correspondent
C a l influence de la temperature et de
C l endommagement (travail inverse de ce qui est
C fait dans le procedure increme).
C================================================
C
********* materiau dependant de la temperature ***********************
         CALL DERTRA(NYOG1,YOG1,T0,YUNG0,YUNGV0,TO,TO)
         CALL DERTRA(NYNU1,YNU1,T0,ENU0,ENUV0,TO,TO)
         XMAT(1)=YUNG0
         XMAT(2)=ENU0
C------------------------------------------------
C Calcul de la matrice de Hooke inverse DD a t=t0
C------------------------------------------------
         CALL ELAST4(2,IFOURB,VAR0,NVARI,XMAT,NCOMAT,YUNGV0,ENUV0,
     &               XCARB,ICARA,MFR1,NSTRS1,DD,DDV,KERRE,2,ITHHER)
C
********* materiau dependant de la temperature ***********************
         CALL DERTRA(NYOG1,YOG1,TF,YUNG1,YUNGV1,TO,TO)
         CALL DERTRA(NYNU1,YNU1,TF,ENU1,ENUV1,TO,TO)
         XMAT(1)=YUNG1
         XMAT(2)=ENU1
C---------------------------------------------------
C Calcul de la matrice de Hooke inverse DDINV a T=TF
C---------------------------------------------------
         CALL ELAST4(2,IFOURB,VAR0,NVARI,XMAT,NCOMAT,YUNGV1,ENUV1,
     &               XCARB,ICARA,MFR1,NSTRS1,DDINV,DDV,KERRE,2,ITHHER)
***********************************************************************
C
         CALL DERTRA(NYALFT1,YALFT1,T0,ALFA0,ALFAV0,TO,TO)
         CALL DERTRA(NYALFT1,YALFT1,TF,ALFAF,ALFAVF,TO,TO)
         CTEPS=ALFA0*(T0-TREF)-ALFAF*(TF-TREF)
C
         CALL ZDANUL(DSPT,NSTRS1)
         AA=1.D0
         DO I=1,NSTRS1
           DSPT(I)=DEPST(I)
           IF (I.GT.3) AA=0.D0
           DO J=1,NSTRS1
             DSPT(I)=DSPT(I)+(DDINV(I,J)*SIG0(J))
             DSPT(I)=DSPT(I)-(DD(I,J)*SIG0(J))
           ENDDO
           DSPT(I)=DSPT(I)-(AA*CTEPS)
           DSPT(I)=DSPT(I)/DT
         ENDDO
       ENDIF
       DO I=1,NVARI
         VAR(I)=VAR0(I)
       ENDDO
       IF (NVARI.LT.100) THEN
          DO I=NVARI+1,100
            VAR(I)=0.0D0
          ENDDO
       ENDIF
 
 
C
C ---------------------------------------------------------------------
       INV = 0
       NSSINC = 0
       nitera = 1
       nopri=0
C ---------------------------------------------------------------------
C       DEBUT DES ITERATIONS EN SSINCREMENTS  /FIN SI DTLEFT = 0
C ---------------------------------------------------------------------
  70   CONTINUE
       iforce=.false.
C       if(dtaumi.gt.TAU) THEN
       if(dtaumi.gt.TAU.and.nssinc.gt.400) THEN
        tau=min(dtaumi,dtleft)
        iforce=.true.
       endif
       NSSINC = NSSINC + 1
       nopri = nopri + nitera - 1
       IF (NSSINC.GT.msoupa) THEN
        DTLIBR=.FALSE.
        GOTO 6543
C       GOTO 999
       ENDIF
C
C---------------------------------------------------------------------
C      START OF CALCULATIONS
C_____________________________________________________________________
       IF (MFR1.EQ.17.AND.INPLAS.NE.19) GOTO 999
       IF (MFR1.NE.3.AND.MFR1.NE.17.AND.INPLAS.EQ.17) THEN
         CALL XXPT1(SIG,EPSV,VAR,EVP1,VARP1,XPM1,XMAT,NSTRS1,NVARI,
     &              NCOMAT,MFR1)
       ELSE IF (MFR1.NE.3.AND.MFR1.NE.17.AND.((INPLAS.GE.19.AND.
     &          INPLAS.LE.24).OR.INPLAS.EQ.61.OR.INPLAS.EQ.107)) THEN
C-----------------------------------------------------------
C Mise a jour eventuelle de la nouvelle temperature TI1 et
C de la nouvelle densite de fissions FII1
C-----------------------------------------------------------
C
         DELTAT=TPOINT*TAU
         TI1=TI0+DELTAT
         DELTAF=FPOINT*TAU
         FII1=FII0+DELTAF
         CALL INCRE1(TAU,SIG,EPSV,VAR,EVP1,VARP1,XMAT,NSTRS1,NVARI,
     &               INPLAS,NCOMAT,MFR1,FII0,FII1,TI1)
       ELSE IF (MFR1.NE.3.AND.MFR1.NE.17.AND.INPLAS.EQ.25) THEN
         CALL INCRE2(SIG,VAR,EVP1,VARP1,XMAT,XMAT0,NSTRS1,MFR1,NVARI,
     &               NCOMAT,VALMAT,VALMA0,DT0)
       ELSE IF ((MFR1.EQ.1.OR.MFR1.EQ.31).AND.INPLAS.EQ.130) THEN
         CALL INCREP(SIG,VAR,EVP1,VARP1,XMAT,NSTRS1,MFR1,
     &               NVARI,NCOMAT)
       ELSE IF ((MFR1.EQ.1.OR.MFR1.EQ.31).AND.
     &(INPLAS.EQ.136.OR.INPLAS.EQ.137.OR.
     &INPLAS.EQ.138.OR.INPLAS.EQ.139)) THEN
         CALL INCVAR(SIG,VAR,EVP1,VARP1,XMAT,NSTRS1,MFR1,
     &               NVARI,NCOMAT,INPLAS)
       ELSE IF (MFR1.NE.3.AND.MFR1.NE.17.AND.INPLAS.EQ.76) THEN
         CALL INCRA2(SIG,VAR,EVP1,VARP1,XMAT,NSTRS1,MFR1,NVARI,NCOMAT)
       ELSE IF (MFR1.NE.3.AND.MFR1.NE.17.AND.INPLAS.EQ.77) THEN
         CALL INCRB2(SIG,VAR,EVP1,VARP1,XMAT,NSTRS1,MFR1,NVARI,NCOMAT)
       ELSE IF (MFR1.NE.3.AND.MFR1.NE.17.AND.INPLAS.EQ.53) THEN
         CALL INCRE5(SIG,VAR,EVP1,VARP1,XMAT,NSTRS1,MFR1,NVARI,NCOMAT)
       ELSE IF (MFR1.NE.3.AND.INPLAS.EQ.63) THEN
         CALL INCRE7(SIG,VAR,DSPT,EVP1,VARP1,XMAT,NSTRS1,MFR1,NVARI,
     &               NCOMAT,IFOURB)
       ELSE IF (MFR1.EQ.3.AND.((INPLAS.GE.19.AND.INPLAS.LE.24)
     &   .OR.INPLAS.EQ.61)) THEN
         CALL INCRE3(TAU,SIG,EPSV,VAR,XMAT,EVP1,VARP1,ALFA,NSTRS1,
     &               NVARI,INPLAS,NCOMAT)
       ELSE IF (MFR1.EQ.3.AND.INPLAS.EQ.25) THEN
         CALL INCRE4(SIG,VAR,EVP1,VARP1,XMAT,ALFA,NSTRS1,NVARI,NCOMAT)
       ELSE IF (MFR1.EQ.3.AND.INPLAS.EQ.76) THEN
         CALL INCRA4(SIG,VAR,EVP1,VARP1,XMAT,ALFA,NSTRS1,NVARI,NCOMAT)
       ELSE IF (MFR1.EQ.3.AND.INPLAS.EQ.77) THEN
         CALL INCRB4(SIG,VAR,EVP1,VARP1,XMAT,ALFA,NSTRS1,NVARI,NCOMAT)
       ELSE IF (MFR1.EQ.3.AND.INPLAS.EQ.53) THEN
         CALL INCRE6(SIG,VAR,EVP1,VARP1,XMAT,ALFA,NSTRS1,NVARI,NCOMAT)
       ELSE IF (MFR1.NE.3.AND.(INPLAS.EQ.85.OR.INPLAS.EQ.86.
     &          OR.INPLAS.EQ.84.OR.INPLAS.EQ.102)) THEN
         CALL DEVFLO(INPLAS,SIG,EPSV,VAR,XMAT,NCOMAT,NSTRS1,NVARI,EVP1,
     &               VARP1,TAU)
       ELSE IF (MFR1.NE.3.AND.INPLAS.EQ.70) THEN
         CALL INCRE8(SIG,VAR,T0,TF,EVP1,VARP1,XMAT,NSTRS1,MFR1,NVARI,
     &               NCOMAT)
C
       ELSE IF (MFR1.NE.3.AND.MFR1.NE.17.AND.INPLAS.EQ.29) THEN
         if (nssinc.eq.1) then
C
C-----------------------------------------------------------
C Initialisation du tableau XMAT(NCOMAT) a l'instant t=t0
C ou le materiau est a la temperature TI0 comprise dans [TINF,TSUP]
C-----------------------------------------------------------
           CALL INITT(YOG,NYOG,YNU,NYNU,YN,NYN,YM,NYM,YKK,NYKK,YALFA1,
     &                NYALF1,YBETA1,NYBET1,YR,NYR,YA,NYA,YALFA,NYALFA,
     &                YRHO,NYRHO,SIGY,NSIGY,YKX,NYKX,NKX,NNKX,XMAT,
     &                NCOMAT,TI0,TINF,TSUP,ITEST,TRUC,NCOURB)
C
           IF (ITHHER.EQ.2) THEN
C********** materiau dependant de la temperature **********************
C---------- Initialisation du tableau XMAT1(NCOMAT) a T=TINF
             CALL INITT(YOG,NYOG,YNU,NYNU,YN,NYN,YM,NYM,YKK,NYKK,
     &                  YALFA1,NYALF1,YBETA1,NYBET1,YR,NYR,YA,NYA,
     &                  YALFA,NYALFA,YRHO,NYRHO,SIGY,NSIGY,YKX,NYKX,
     &                  NKX,NNKX,XMAT1,NCOMAT,TINF,TO,TO,ITEST,TRUC,
     &                  NCOURB)
C---------- Initialisation du tableau XMAT2(NCOMAT) a T=TSUP
             CALL INITT(YOG,NYOG,YNU,NYNU,YN,NYN,YM,NYM,YKK,NYKK,
     &                  YALFA1,NYALF1,YBETA1,NYBET1,YR,NYR,YA,NYA,
     &                  YALFA,NYALFA,YRHO,NYRHO,SIGY,NSIGY,YKX,NYKX,
     &                  NKX,NNKX,XMAT2,NCOMAT,TSUP,TO,TO,ITEST,TRUC,
     &                  NCOURB)
C**********************************************************************
           ENDIF
C--------------------------------------------------------------
C Calcul de la derivee de la dilatation thermique /temps a t=t0
C--------------------------------------------------------------
           CALL DERTRA(NYALFA,YALFA,TI0,ALFA0,ALFAV0,TO,TO)
           CALL ZDANUL(EPSTHD,NSTRS1)
           CTH=(ALFAV0*TPOINT*(TI0-TREF))+(ALFA0*TPOINT)
           DO I=1,3
             EPSTHD(I)=CTH
           ENDDO
C
           CALL XXCREE(TAU,SIG,EPSV,VAR,SIGP1,EVP1,VARP1,EPSTHD,DSPT,
     &                 XMAT,XMAT1,XMAT2,YKX,NKX,NSTRS1,NVARI,NCOMAT,
     &                 NYKX,NNKX,DD,DDV,DDINV,YOG,NYOG,YNU,NYNU,MFR1,
     &                 XCARB,ICARA,IFOURB,2,TI0,TPOINT,TINF,TSUP,ITEST,
     &                 ITHHER,TRUC,NCOURB)
c         DO i= 1,6
c           WRITE(6,*) 'sigp1(',i,')= ', sigp1(i)
c         ENDDO
c         DO i= 1,6
c           WRITE(6,*) 'evp1(',i,')= ', evp1(i)
c         ENDDO
c         DO i= 1,6
c           WRITE(6,*) 'varp1(',i,')= ', varp1(i)
c         ENDDO
           CALL ESTITO(SIG,NSTRS1,VAR,NVARI,YKX,NYKX,NKX,NNKX,XMAT,
     &                 NCOMAT,TI0,TAUX,TRUC,NCOURB)
         else
           do i=1,nstrs
             sigp1(i)=sigp4(i)
             evp1(i)=evp4(i)
           enddo
           do i=1,nvari
             varp1(i)=varp4(i)
           enddo
         endif
C
       ELSE IF (MFR1.NE.3.AND.MFR1.NE.17.AND.INPLAS.EQ.142) THEN
         if (nssinc.eq.1) then
C
C-----------------------------------------------------------
C Initialisation du tableau XMAT(NCOMAT) a l'instant t=t0
C ou le materiau est a la temperature TI0 comprise dans [TINF,TSUP]
C-----------------------------------------------------------
           CALL INITT1(YOG1,NYOG1,YNU1,NYNU1,YN1,NYN1,YM1,NYM1,YKK1,
     &                NYKK1,YALF2,NYALF2,YBET2,NYBET2,YR1,NYR1,YA1,NYA1,
     &                YQ1,NYQ1,YALFT1,NYALFT1,YRHO1,NYRHO1,SIGY1,NSIGY1,
     &                XMAT,NCOMAT,TI0,TINF,TSUP,ITEST)
C
           IF (ITHHER.EQ.2) THEN
C********** materiau dependant de la temperature **********************
C---------- Initialisation du tableau XMAT1(NCOMAT) a T=TINF
             CALL INITT1(YOG1,NYOG1,YNU1,NYNU1,YN1,NYN1,YM1,NYM1,YKK1,
     &                  NYKK1,YALF2,NYALF2,YBET2,NYBET2,YR1,NYR1,YA1,
     &                  NYA1,YQ1,NYQ1,YALFT1,NYALFT1,YRHO1,NYRHO1,SIGY1,
     &                  NSIGY1,XMAT1,NCOMAT,TINF,TO,TO,ITEST)
C---------- Initialisation du tableau XMAT2(NCOMAT) a T=TSUP
             CALL INITT1(YOG1,NYOG1,YNU1,NYNU1,YN1,NYN1,YM1,NYM1,YKK1,
     &                  NYKK1,YALF2,NYALF2,YBET2,NYBET2,YR1,NYR1,YA1,
     &                  NYA1,YQ1,NYQ1,YALFT1,NYALFT1,YRHO1,NYRHO1,SIGY1,
     &                  NSIGY1,XMAT1,NCOMAT,TSUP,TO,TO,ITEST)
C**********************************************************************
           ENDIF
C--------------------------------------------------------------
C Calcul de la derivee de la dilatation thermique /temps a t=t0
C--------------------------------------------------------------
           CALL DERTRA(NYALFT1,YALFT1,TI0,ALFA0,ALFAV0,TO,TO)
           CALL ZDANUL(EPSTHD,NSTRS1)
           CTH=(ALFAV0*TPOINT*(TI0-TREF))+(ALFA0*TPOINT)
           DO I=1,3
             EPSTHD(I)=CTH
           ENDDO
            CALL XXCRE1(TAU,SIG,EPSV,VAR,SIGP1,EVP1,VARP1,EPSTHD,DSPT,
     &                 XMAT,XMAT1,XMAT2,NSTRS1,NVARI,NCOMAT,DD,DDV,
     &                 DDINV,DDINVp,YOG1,NYOG1,YNU1,NYNU1,MFR1,XCARB,
     &                 ICARA,IFOURB,2,TI0,TPOINT,TINF,TSUP,ITEST,
     &                 ITHHER,VART,IB,IGAU,kerre)
c         DO i= 1,6
c           WRITE(6,*) 'sigp1(',i,')= ', sigp1(i)
c         ENDDO
c         DO i= 1,6
c           WRITE(6,*) 'evp1(',i,')= ', evp1(i)
c         ENDDO
c         DO i= 1,6
c           WRITE(6,*) 'varp1(',i,')= ', varp1(i)
c         ENDDO
           VAR(8)=VART(8)
           CALL ESTIT1(SIG,NSTRS1,VAR,NVARI,XMAT,NCOMAT,TI0,TAUX)
         else
           do i=1,nstrs
             sigp1(i)=sigp4(i)
             evp1(i)=evp4(i)
           enddo
           do i=1,nvari
             varp1(i)=varp4(i)
           enddo
         endif
       ELSE IF (MFR1.EQ.17.AND.INPLAS.EQ.19) THEN
         CALL TUFINC(TAU,SIG,EPSV,VAR,XMAT,XCARB,EVP1,VARP1,NSTRS1,
     &               NVARI,INPLAS,NCOMAT,KERREU1)
       ELSE IF ((MFR1.EQ.1.OR.MFR1.EQ.31.OR.MFR1.EQ.33)
     &.AND.INPLAS.EQ.44) THEN
         CALL POUDRA(SIG,EPSV,VAR,EVP1,VARP1,XMAT,NSTRS1,NVARI,NCOMAT,
     &               KERRE)
       ELSE IF ((MFR1.EQ.1.OR.MFR1.EQ.31.OR.MFR1.EQ.33)
     &.AND.INPLAS.EQ.45) THEN
         CALL POUDRB(SIG,EPSV,VAR,EVP1,VARP1,XMAT,NSTRS1,NVARI,NCOMAT,
     &               KERRE)
C------
       ELSE IF ((MFR1.EQ.1.OR.MFR1.EQ.31.OR.MFR1.EQ.33)
     &          .AND.INPLAS.EQ.165) THEN
C        Chaboche
         DO I=1,NCOMAT
            ZMAT(I) = XMAT(I)
         ENDDO
         CALL INCRE9(SIG,VAR,ZMAT,NSTRS1,EVP1,VARP1)
C------
       ENDIF
C
       NITERA  = 0
C --------------------------------------------------------------------
C      DEBUT DES ITERATIONS SUR TAU OPTIMAL /FIN SI RA PETIT
C --------------------------------------------------------------------
  80   CONTINUE
       iforce=.false.
C       WRITE(6,*) 'NITERA', nitera
*        if( tau.lt.dtaumi) then
        if( tau.lt.dtaumi.and.nssinc.gt.400) then
          tau=min(dtaumi,dtleft)
          iforce=.true.
        endif
       IF (MFR1.EQ.3) GOTO 150
         IF (MFR1.EQ.17) GOTO 210
         IF (INPLAS.EQ.17) THEN
           tau2=tau*0.5d0
           CALL AVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,EVP1,VARP1,
     &                XPM1,DSPT,XMAT,NSTRS1,NVARI,NCOMAT,IFOURB,IVTEST,
     &                MFR1)
           CALL XXPT1(SIG1,EPSV1,VAR1,EVP2,VARP2,XPM2,XMAT,NSTRS1,
     &                NVARI,NCOMAT,MFR1)
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
             XPM2(I) = 0.5D0*(XPM1(I)+XPM2(I))
           ENDDO
           DO I=1,4+NSTRS1
             VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
           ENDDO
           CALL AVAN1(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,EVP2,VARP2,
     &                XPM2,DSPT,XMAT,NSTRS1,NVARI,NCOMAT,IFOURB,IVTEST,
     &                MFR1)
           CALL XXPT1(SIG12,EPSV12,VAR12,EVP3,VARP3,XPM3,XMAT,NSTRS1,
     &                NVARI,NCOMAT,MFR1)
           CALL AVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,EVP3,
     &                VARP3,XPM3,DSPT,XMAT,NSTRS1,NVARI,NCOMAT,IFOURB,
     &                IVTEST,MFR1)
           CALL XXPT1(SIG13,EPSV13,VAR13,EVP4,VARP4,XPM4,XMAT,NSTRS1,
     &                NVARI,NCOMAT,MFR1)
           DO I=1,NSTRS1
             EVP4(I) = 0.5D0*(EVP3(I)+EVP4(I))
             XPM4(I) = 0.5D0*(XPM3(I)+XPM4(I))
           enddo
           DO I=1,4+NSTRS1
             VARP4(I) = 0.5D0*(VARP3(I)+VARP4(I))
           enddo
           CALL AVAN1(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,EVP4,
     &                VARP4,XPM4,DSPT,XMAT,NSTRS1,NVARI,NCOMAT,IFOURB,
     &                IVTEST,MFR1)
           CALL XXPT1(SIGf,EPinf,VARf,EVP4,VARP4,XPM4,XMAT,NSTRS1,
     &                NVARI,NCOMAT,MFR1)
           DO I=1,NSTRS1
             EVP2(I) = (EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
             XPM2(I) = (XPM1(I)+XPM4(I))/6.d0+XPM3(I)*2.d0/3.d0
           enddo
           DO I=1,4+NSTRS1
             VARP2(I) = (VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
           enddo
           CALL AVAN1(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,EVP2,VARP2,XPM2,
     &                DSPT,XMAT,NSTRS1,NVARI,NCOMAT,IFOURB,IVTEST,MFR1)
C---------
         ELSE IF (INPLAS.EQ.44) THEN
           CALL AVANP(TAU,SIG,EPSV,VAR,SIGf,EPinf,VARf,DSPT,EVP1,
     &                VARP1,XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
           CALL POUDRA(SIGf,EPinf,VARf,EVP2,VARP2,XMAT,NSTRS1,NVARI,
     &                 NCOMAT,KERRE)
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
           ENDDO
           VARP2(1) = 0.5D0*(VARP1(1)+VARP2(1))
           CALL AVANP(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                VARP2,XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
*           tau2=tau*0.5d0
*           CALL AVANP(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
*     &                VARP1,XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
*           CALL POUDRA(SIG1,EPSV1,VAR1,EVP2,VARP2,XMAT,NSTRS1,NVARI,
*     &                 NCOMAT,KERRE)
*           DO I=1,NSTRS1
*             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
*           ENDDO
*           VARP2(1) = 0.5D0*(VARP1(1)+VARP2(1))
*           CALL AVANP(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
*     &                VARP2,XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
*          CALL POUDRA(SIG12,EPSV12,VAR12,EVP3,VARP3,XMAT,NSTRS1,NVARI,
*     &                 NCOMAT,KERRE)
*           CALL AVANP(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
*     &                EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
*           CALL POUDRA(SIG13,EPSV13,VAR13,EVP4,VARP4,XMAT,NSTRS1,NVARI,
*     &                 NCOMAT,KERRE)
*           DO I=1,NSTRS1
*             EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
*           enddo
*           VARP4(1) = 0.5d0*(VARP3(1)+VARP4(1))
*           CALL AVANP(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
*     &                EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
*           CALL POUDRA(SIGf,EPinf,VARf,EVP4,VARP4,XMAT,NSTRS1,NVARI,
*     &                 NCOMAT,KERRE)
*           DO I=1,NSTRS1
*             EVP2(I) = (EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
*           enddo
*           VARP2(1) = (VARP1(1)+VARP4(1))/6.d0+VARP3(1)*2.d0/3.d0
*           CALL AVANP(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,VARP2,
*     &                XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
*C---------
         ELSE IF (INPLAS.EQ.45) THEN
           CALL AVANP(TAU,SIG,EPSV,VAR,SIGf,EPinf,VARf,DSPT,EVP1,
     &                VARP1,XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
           CALL POUDRB(SIGf,EPinf,VARf,EVP2,VARP2,XMAT,NSTRS1,NVARI,
     &                 NCOMAT,KERRE)
           DO I=1,NSTRS1
             EVP2(I) = 0.5d0*(EVP1(I)+EVP2(I))
           ENDDO
           VARP2(1) = 0.5d0*(VARP1(1)+VARP2(1))
           CALL AVANP(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                VARP2,XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
C---------
C    CALCUL DE LA TAILLE DE GRAIN
C---------
           CALL GRAIN(TAU,EVP1,EVP2,SIG,SIG1,VAR,VAR1,XMAT,NSTRS1,
     &                NVARI,KERRE)
 
 
*           tau2=tau*0.5d0
*           CALL AVANP(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
*     &                VARP1,XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
*           CALL POUDRB(SIG1,EPSV1,VAR1,EVP2,VARP2,XMAT,NSTRS1,NVARI,
*     &                 NCOMAT,KERRE)
*           DO I=1,NSTRS1
*             EVP2(I) = 0.5d0*(EVP1(I)+EVP2(I))
*           ENDDO
*           VARP2(1) = 0.5d0*(VARP1(1)+VARP2(1))
*           CALL AVANP(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
*     &                VARP2,XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
*C---------
*C    CALCUL DE LA TAILLE DE GRAIN
*C---------
*           CALL GRAIN(TAU2,EVP1,EVP2,SIG,SIG12,VAR,VAR12,XMAT,NSTRS1,
*     &                NVARI,KERRE)
*c
*           CALL POUDRB(SIG12,EPSV12,VAR12,EVP3,VARP3,XMAT,NSTRS1,NVARI,
*     &                 NCOMAT,KERRE)
*           CALL AVANP(TAU2,SIG12,EPSV12,VAR12,SIG13,EPSV13,VAR13,DSPT,
*     &                EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
*           CALL POUDRB(SIG13,EPSV13,VAR13,EVP4,VARP4,XMAT,NSTRS1,NVARI,
*     &                 NCOMAT,KERRE)
*           DO I=1,NSTRS1
*             EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
*           enddo
*           VARP4(1) = 0.5d0*(VARP1(1)+VARP2(1))
*           CALL AVANP(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
*     &                EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
*           CALL POUDRB(SIGf,EPinf,VARf,EVP4,VARP4,XMAT,NSTRS1,NVARI,
*     &                 NCOMAT,KERRE)
*           DO I=1,NSTRS1
*             EVP2(I) =(EVP3(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
*           enddo
*           VARP2(1) =(VARP1(1)+VARP2(1))/6.d0+VARP3(1)*2.d0/3.d0
*           CALL AVANP(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,VARP2,
*     &                XMAT,NSTRS1,NVARI,IFOURB,NCOMAT)
*C---------
*C    CALCUL DE LA TAILLE DE GRAIN
*C---------
*           CALL GRAIN(TAU2,EVP3,EVP4,SIG12,SIGf,VAR12,VARf,XMAT,NSTRS1,
*     &                NVARI,KERRE)
*C---------
         ELSE IF ((INPLAS.GE.19.AND.INPLAS.LE.24).OR.INPLAS.EQ.61.OR.
     &            INPLAS.EQ.107) THEN
           tau2=tau*0.5d0
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C-----------------------------------------------------------
C Mise a jour eventuelle de la nouvelle temperature TI1 et
C de la nouvelle densite de fissions FII1
C-----------------------------------------------------------
C
           DELTAT=TPOINT*TAU
           TI1=TI0+DELTAT
           DELTAF=FPOINT*TAU
           FII1=FII0+DELTAF
           IF (INPLAS.EQ.107) VAR1(3)=VAR(3)
           CALL INCRE1(TAU,SIG1,EPSV1,VAR1,EVP2,VARP2,XMAT,NSTRS1,
     &                 NVARI,INPLAS,NCOMAT,MFR1,FII0,FII1,TI1)
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
           ENDDO
           IF (INPLAS.EQ.24.OR.INPLAS.EQ.107) THEN
             IF (INPLAS.EQ.107) VARPBU=VARP2(3)
             DO I=1,NVARI
               VARP2(I)=0.5D0*(VARP1(I)+VARP2(I))
             ENDDO
           ELSE
             DO I=1,2*NSTRS1+2
               VARP2(I)= 0.5D0*(VARP1(I) + VARP2(I))
             ENDDO
             DO I=2*NSTRS1+4,NVARI
               VARP2(I)= 0.5D0*(VARP1(I) + VARP2(I))
             ENDDO
           ENDIF
C-----------------------------------------------------------
C Mise a jour eventuelle de la nouvelle temperature TI1 et
C de la nouvelle densite de fissions FII1
C-----------------------------------------------------------
C
           DELTAT=TPOINT*TAU2
           TI12=TI0+DELTAT
           DELTAF=FPOINT*TAU2
           FII12=FII0+DELTAF
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           IF (INPLAS.EQ.107) VAR12(3)=VAR(3)
           CALL INCRE1(TAU2,SIG12,EPSV12,VAR12,EVP3,VARP3,XMAT,NSTRS1,
     &                 NVARI,INPLAS,NCOMAT,MFR1,FII0,FII12,TI12)
C-----------------------------------------------------------
C Mise a jour eventuelle de la nouvelle temperature TI1 et
C de la nouvelle densite de fissions FII1
C-----------------------------------------------------------
C
           TI1=TI12+DELTAT
           FII1=FII12+DELTAF
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE1(TAU2,SIG13,EPSV13,VAR13,EVP4,VARP4,XMAT,NSTRS1,
     &                 NVARI,INPLAS,NCOMAT,MFR1,FII12,FII1,TI1)
           DO I=1,NSTRS1
             EVP4(I) = 0.5D0*(EVP3(I)+EVP4(I))
           enddo
           IF (INPLAS.EQ.24.OR.INPLAS.EQ.107) THEN
             DO I=1,NVARI
               VARP4(I)=0.5D0*(VARP3(I)+VARP4(I))
             enddo
           ELSE
             DO I=1,2*NSTRS1+2
               VARP4(I)= 0.5D0*(VARP3(I) + VARP4(I))
             enddo
             DO  I=2*NSTRS1+4,NVARI
               VARP4(I)= 0.5D0*(VARP3(I) + VARP4(I))
             enddo
           ENDIF
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C-----------------------------------------------------------
C Mise a jour eventuelle de la nouvelle temperature TI1 et
C de la nouvelle densite de fissions FII1
C-----------------------------------------------------------
C
           DELTAT=TPOINT*TAU2
           TI12=TI0+DELTAT
           DELTAF=FPOINT*TAU2
           FII12=FII0+DELTAF
           IF (INPLAS.EQ.107) VARf(3)=VAR(3)
           CALL INCRE1(TAU2,SIGf,EPinf,VARf,EVP4,VARP4,XMAT,NSTRS1,
     &                 NVARI,INPLAS,NCOMAT,MFR1,FII0,FII12,TI12)
           DO I=1,NSTRS1
             EVP2(I) = (EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           enddo
           IF (INPLAS.EQ.24.OR.INPLAS.EQ.107) THEN
             DO I=1,NVARI
               VARP2(I)=(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
             enddo
             IF (INPLAS.EQ.107) VARP2(3)=VARPBU
           ELSE
             DO I=1,2*NSTRS1+2
               VARP2(I)=(VARP1(I) + VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
             enddo
             DO  I=2*NSTRS1+4,NVARI
               VARP2(I)= (VARP1(I) + VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
             enddo
           ENDIF
           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAr1,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           IF(INPLAS.EQ.107) VARf(3)=VAR1(3)
C---------
         ELSE IF (INPLAS.EQ.25) THEN
           tau2=tau*0.5d0
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE2(SIG1,VAR1,EVP2,VARP2,XMAT,XMAT0,NSTRS1,MFR1,
     &                 NVARI,NCOMAT,VALMAT,VALMA0,DT0)
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
           ENDDO
           DO I=1,NVARI
            VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
           ENDDO
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE2(SIG12,VAR12,EVP3,VARP3,XMAT,XMAT0,NSTRS1,MFR1,
     &                 NVARI,NCOMAT,VALMAT,VALMA0,DT0)
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE2(SIG13,VAR13,EVP4,VARP4,XMAT,XMAT0,NSTRS1,MFR1,
     &                 NVARI,NCOMAT,VALMAT,VALMA0,DT0)
           DO I=1,NSTRS1
             EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
           enddo
           DO I=1,NVARI
             VARP4(I) = 0.5d0*(VARP3(I)+VARP4(I))
           enddo
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE2(SIGf,VARf,EVP4,VARP4,XMAT,XMAT0,NSTRS1,MFR1,
     &                 NVARI,NCOMAT,VALMAT,VALMA0,DT0)
           DO I=1,NSTRS1
             EVP2(I) =(EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           enddo
           DO I=1,NVARI
             VARP2(I) =(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
           enddo
           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C---------
         ELSE IF (INPLAS.EQ.130) THEN
           tau2=tau*0.5d0
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCREP(SIG1,VAR1,EVP2,VARP2,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
           ENDDO
           DO I=1,NVARI
            VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
           ENDDO
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCREP(SIG12,VAR12,EVP3,VARP3,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT)
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCREP(SIG13,VAR13,EVP4,VARP4,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
           enddo
           DO I=1,NVARI
             VARP4(I) = 0.5d0*(VARP3(I)+VARP4(I))
           enddo
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCREP(SIGf,VARf,EVP4,VARP4,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP2(I) =(EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           enddo
           DO I=1,NVARI
             VARP2(I) =(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
           enddo
           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C---------
C--------- LOI CHAB_SINH_R, CHAB_SINH_X (Runge Kuta 2-3)
         ELSE IF (INPLAS.EQ.136.OR.INPLAS.EQ.137.OR.
     a            INPLAS.EQ.138.OR.INPLAS.EQ.139) THEN
           tau2=tau*0.5d0
C        write(6,*)'première entrée dans advac'
           CALL ADVAR(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &             DEPST,VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C        write(6,*)'sortie de advac'
C        write(6,*)'seconde entrée dans increp'
           CALL INCVAR(SIG1,VAR1,EVP2,VARP2,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT,INPLAS)
C        write(6,*)'sortie de increp'
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
           ENDDO
           DO I=1,2
            VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
           ENDDO
C        write(6,*)'seconde entrée dans advac'
           CALL ADVAR(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
     &              DEPST,VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C        write(6,*)'sortie de advac'
C        write(6,*)'troisième entrée entrée dans increp'
           CALL INCVAR(SIG12,VAR12,EVP3,VARP3,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT,INPLAS)
C        write(6,*)'sortie de increp'
C        write(6,*)'troisième entrée dans advac'
           CALL ADVAR(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &          EVP3,DEPST,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C        write(6,*)'sortie de advac'
C        write(6,*)'quatrième entrée dans increp'
           CALL INCVAR(SIG13,VAR13,EVP4,VARP4,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT,INPLAS)
C        write(6,*)'sortie de increp'
           DO I=1,NSTRS1
             EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
           enddo
           DO I=1,2
             VARP4(I) = 0.5d0*(VARP3(I)+VARP4(I))
           enddo
C        write(6,*)'quatrième entrée dans advac'
           CALL ADVAR(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
     &          EVP4,DEPST,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C        write(6,*)'sortie de advac'
C        write(6,*)'cinquième entrée dans increp'
           CALL INCVAR(SIGf,VARf,EVP4,VARP4,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT,INPLAS)
C        write(6,*)'sortie de increp'
           DO I=1,NSTRS1
             EVP2(I) =(EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           enddo
           DO I=1,2
             VARP2(I) =(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
           enddo
C        write(6,*)'cinquième entrée dans advac'
           CALL ADVAR(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &          DEPST,VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C        write(6,*)'sortie de advac'
C---------
C---------
         ELSE IF (INPLAS.EQ.76) THEN
           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIGf,EPinf,VARf,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRA2(SIGf,VARf,EVP2,VARP2,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP2(I) = 0.5d0*(EVP1(I)+EVP2(I))
           ENDDO
           DO I=1,NVARI
             VARP2(I) = 0.5d0*(VARP1(I)+VARP2(I))
           ENDDO
           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
*           tau2=tau*0.5d0
*           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
*     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
*           CALL INCRA2(SIG1,VAR1,EVP2,VARP2,XMAT,NSTRS1,MFR1,NVARI,
*     &                 NCOMAT)
*           DO I=1,NSTRS1
*             EVP2(I) = 0.5d0*(EVP1(I)+EVP2(I))
*           ENDDO
*           DO I=1,NVARI
*             VARP2(I) = 0.5d0*(VARP1(I)+VARP2(I))
*           ENDDO
*           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
*     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
*           CALL INCRA2(SIG12,VAR12,EVP3,VARP3,XMAT,NSTRS1,MFR1,NVARI,
*     &                 NCOMAT)
*           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
*     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
*           CALL INCRA2(SIG13,VAR13,EVP4,VARP4,XMAT,NSTRS1,MFR1,NVARI,
*     &                 NCOMAT)
*           DO I=1,NSTRS1
*             EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
*           enddo
*           DO I=1,NVARI
*             VARP4(I) = 0.5d0*(VARP3(I)+VARP4(I))
*           enddo
*           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
*     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
*           CALL INCRA2(SIGf,VARf,EVP4,VARP4,XMAT,NSTRS1,MFR1,NVARI,
*     &                 NCOMAT)
*           DO I=1,NSTRS1
*             EVP2(I) =(EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
*           enddo
*           DO I=1,NVARI
*             VARP2(I) = (VARP3(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
*           enddo
*           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
*     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C---------
         ELSE IF (INPLAS.EQ.77) THEN
           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIGf,EPinf,VARf,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRB2(SIGf,VARf,EVP2,VARP2,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP2(I) = 0.5d0*(EVP1(I)+EVP2(I))
           enddo
           DO I=1,NVARI
             VARP2(I) = 0.5d0*(VARP1(I)+VARP2(I))
           enddo
           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
 
*           tau2=tau*0.5d0
*           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
*     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
*           CALL INCRB2(SIG1,VAR1,EVP2,VARP2,XMAT,NSTRS1,MFR1,NVARI,
*     &                 NCOMAT)
*           DO I=1,NSTRS1
*             EVP2(I) = 0.5d0*(EVP1(I)+EVP2(I))
*           enddo
*           DO I=1,NVARI
*             VARP2(I) = 0.5d0*(VARP1(I)+VARP2(I))
*           enddo
*           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
*     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
*           CALL INCRB2(SIG12,VAR12,EVP3,VARP3,XMAT,NSTRS1,MFR1,NVARI,
*     &                 NCOMAT)
*           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
*     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
*           CALL INCRB2(SIG13,VAR13,EVP4,VARP4,XMAT,NSTRS1,MFR1,NVARI,
*     &                 NCOMAT)
*           DO I=1,NSTRS1
*             EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
*           enddo
*           DO I=1,NVARI
*             VARP4(I) = 0.5d0*(VARP3(I)+VARP4(I))
*           enddo
*           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
*     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
*           CALL INCRB2(SIGf,VARf,EVP4,VARP4,XMAT,NSTRS1,MFR1,NVARI,
*     &                 NCOMAT)
*           DO I=1,NSTRS1
*             EVP2(I) = (EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
*           enddo
*           DO I=1,NVARI
*             VARP2(I) =(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
*           enddo
*           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
*     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C---------
         ELSE IF (INPLAS.EQ.53) THEN
           tau2=tau*0.5d0
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE5(SIG1,VAR1,EVP2,VARP2,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
           ENDDO
           DO I=1,NVARI
             VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
           ENDDO
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE5(SIG12,VAR12,EVP3,VARP3,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT)
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE5(SIG13,VAR13,EVP4,VARP4,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
           enddo
           DO I=1,NVARI
             VARP4(I) = 0.5d0*(VARP3(I)+VARP4(I))
           enddo
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE5(SIGf,VARf,EVP4,VARP4,XMAT,NSTRS1,MFR1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
              EVP2(I) =(EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           enddo
           DO I=1,NVARI
             VARP2(I) =(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
           enddo
           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C---------
         ELSE IF (INPLAS.EQ.85.OR.INPLAS.EQ.86.OR.INPLAS.EQ.84.OR.
     &       INPLAS.EQ.102) THEN
           tau2=tau*0.5d0
           CALL ADVFLO(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT,
     &                 MFR1)
           CALL DEVFLO(INPLAS,SIG1,EPSV1,VAR1,XMAT,NCOMAT,NSTRS1,NVARI,
     &                 EVP2,VARP2,TAU2)
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
           ENDDO
           DO I=1,NVARI
             VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
           ENDDO
C       >>>>  SI JE SUIS SUR LE PREMIER PAS &lt;&lt;&lt;&lt;
C       et que j'ai calcule le resultat a la main
C       _________________________________________
           IF (VAR(12).GT.1.D0) THEN
C
             DO I=1,NVARI
               VARP2(I) = VARP1(I)
             ENDDO
             DO I=1,NSTRS1
               EVP2(I) = EVP1(I)
             ENDDO
           ENDIF
C
           CALL ADVFLO(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT,
     &                 MFR1)
           CALL DEVFLO(INPLAS,SIG12,EPSV12,VAR12,XMAT,NCOMAT,NSTRS1,
     &                 NVARI,EVP3,VARP3,TAU2)
           CALL ADVFLO(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
     &                 NCOMAT,MFR1)
           CALL DEVFLO(INPLAS,SIG13,EPSV13,VAR13,XMAT,NCOMAT,NSTRS1,
     &                 NVARI,EVP4,VARP4,TAU2)
C           print*,'==>4',VARP4(2),VARP4(3)
           DO I=1,NSTRS1
              EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
           enddo
C
           DO I=1,NVARI
              VARP4(I) = 0.5d0*(VARP3(I)+VARP4(I))
           enddo
           CALL ADVFLO(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
     &                 NCOMAT,MFR1)
           CALL DEVFLO(INPLAS,SIGf,EPinf,VARf,XMAT,NCOMAT,NSTRS1,NVARI,
     &                 EVP4,VARP4,TAU2)
           DO I=1,NSTRS1
              EVP2(I) = (EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           enddo
C
           DO I=1,NVARI
              VARP2(I) = (VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
           enddo
C            print*,'==>5',VARP2(2),VARP2(3)
           CALL ADVFLO(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT,
     &                 MFR1)
C----------
         ELSE IF (INPLAS.EQ.63) THEN
           tau2=tau*0.5d0
           CALL ADVDDI(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT,
     &                 MFR1)
           CALL INCRE7(SIG1,VAR1,DSPT,EVP2,VARP2,XMAT,NSTRS1,MFR1,
     &                 NVARI,NCOMAT,IFOURB)
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
           ENDDO
           DO I=1,NVARI
             VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
           ENDDO
           CALL ADVDDI(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT,
     &                 MFR1)
           CALL INCRE7(SIG12,VAR12,DSPT,EVP3,VARP3,XMAT,NSTRS1,MFR1,
     &                 NVARI,NCOMAT,IFOURB)
           CALL ADVDDI(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
     &                 NCOMAT,MFR1)
           CALL INCRE7(SIG13,VAR13,DSPT,EVP4,VARP4,XMAT,NSTRS1,MFR1,
     &                 NVARI,NCOMAT,IFOURB)
           DO I=1,NSTRS1
              EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
           enddo
           DO I=1,NVARI
              VARP4(I) = 0.5d0*(VARP3(I)+VARP4(I))
           enddo
           CALL ADVDDI(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
     &                 NCOMAT,MFR1)
           CALL INCRE7(SIGf,VARf,DSPT,EVP4,VARP4,XMAT,NSTRS1,MFR1,
     &                 NVARI,NCOMAT,IFOURB)
           DO I=1,NSTRS1
              EVP2(I) =(EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           enddo
           DO I=1,NVARI
              VARP2(I) =(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
           enddo
           CALL ADVDDI(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT,
     &                 MFR1)
C--------
         ELSE IF (INPLAS.EQ.70) THEN
           tau2=tau*0.5d0
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           t12=(t0+tf)*0.5d0
           CALL INCRE8(SIG1,VAR1,T0,T12,EVP2,VARP2,XMAT,NSTRS1,MFR1,
     &                 NVARI,NCOMAT)
C
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
           ENDDO
           DO I=1,NVARI
             VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
           ENDDO
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE8(SIG12,VAR12,T12,TF,EVP3,VARP3,XMAT,NSTRS1,MFR1,
     &                 NVARI,NCOMAT)
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE8(SIG13,VAR13,T12,TF,EVP4,VARP4,XMAT,NSTRS1,MFR1,
     &                 NVARI,NCOMAT)
           DO I=1,NSTRS1
             EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
           enddo
           DO I=1,NVARI
             VARP4(I) = 0.5d0*(VARP3(I)+VARP4(I))
           enddo
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VArf,DSPT,
     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE8(SIGf,VARf,T12,TF,EVP4,VARP4,XMAT,NSTRS1,MFR1,
     &                 NVARI,NCOMAT)
           DO I=1,NSTRS1
             EVP2(I) =(EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           enddo
           DO I=1,NVARI
             VARP2(I) =(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
           enddo
           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAr1,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C---------
         ELSE IF (INPLAS.EQ.29) THEN
C
 143       TAU2=0.5D0*TAU
C           write(6,*) 'tau=',tau
C           write(6,*) 'NSSINC=',NSSINC
*      if (ib.eq.1.and.igau.eq.1)write (6,*) ' tau '  ,tau
           CALL AVANXX(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,SIGP1,EVP1,
     &                 VARP1,NSTRS1,NVARI)
c         write(6,*) 'varp1(3)=',varp1(3)
c       Do i=1,6
c         write(6,*) 'sig1(',i,')=',sig1(i)
c       enddo
c       Do i=1,6
c         write(6,*) 'epsv1(',i,')=',epsv1(i)
c       enddo
c       Do i=1,7
c         write(6,*) 'var1(',i,')=',var1(i)
c       enddo
*       if (ib.eq.1.and.igau.eq.1)
*     $ write(6,*)'sig1(1) sig2 sig3',sig1(1),sig1(2),sig1(3)
           aap = MAX(ABS(SIG1(1)-SIG(1)),ABS(SIG1(2)-SIG(2)))
           aap = max ( aap,ABS (Sig1(3)-SIG(3)))
           IF ( aap . gt . XMAX * 5.) THEN
*       write(6,*)'sig1(1) sig2 sig3',sig1(1),sig1(2),sig1(3)
             rap = aap / xmax
             TAU= TAU / rap
*             if (ib.eq.1.and.igau.eq.1)write(6,*)'rap tau'
             go  to 143
*             do I=1,nstrs
*             sig1(i)=XMAX*100.
*             sigf(I)=xmax*200.
*             enddo
*             go to 250
           endif
C
C-----------------------------------------------------------
C Mise a jour eventuelle de la nouvelle temperature TI1
C-----------------------------------------------------------
C
           DELTAT=TPOINT*TAU2
           TI12=TI0+DELTAT
C Initialisation du tableau XMAT(NCOMAT) a l'instant t=t+TAU/2
C ou le materiau est a la temperature TI12 comprise dans [TINF,TSUP]
C-----------------------------------------------------------
           CALL INITT(YOG,NYOG,YNU,NYNU,YN,NYN,YM,NYM,YKK,NYKK,YALFA1,
     &                NYALF1,YBETA1,NYBET1,YR,NYR,YA,NYA,YALFA,NYALFA,
     &                YRHO,NYRHO,SIGY,NSIGY,YKX,NYKX,NKX,NNKX,XMAT,
     &                NCOMAT,TI12,TINF,TSUP,ITEST,TRUC,NCOURB)
C
           IF (ITHHER.EQ.2) THEN
C********** materiau dependant de la temperature **********************
C---------- Initialisation du tableauXMAT1F(NCOMAT) a T=TINF
             CALL INITT(YOG,NYOG,YNU,NYNU,YN,NYN,YM,NYM,YKK,NYKK,
     &                  YALFA1,NYALF1,YBETA1,NYBET1,YR,NYR,YA,NYA,
     &                  YALFA,NYALFA,YRHO,NYRHO,SIGY,NSIGY,YKX,NYKX,
     &                  NKX,NNKX,XMAT1,NCOMAT,TINF,TO,TO,ITEST,TRUC,
     &                  NCOURB)
C---------- Initialisation du tableauXMAT2P(NCOMAT) a T=TSUP
             CALL INITT(YOG,NYOG,YNU,NYNU,YN,NYN,YM,NYM,YKK,NYKK,
     &                  YALFA1,NYALF1,YBETA1,NYBET1,YR,NYR,YA,NYA,
     &                  YALFA,NYALFA,YRHO,NYRHO,SIGY,NSIGY,YKX,NYKX,
     &                  NKX,NNKX,XMAT2,NCOMAT,TSUP,TO,TO,ITEST,TRUC,
     &                  NCOURB)
C**********************************************************************
           ENDIF
C------------------------------------------------------------------
C Calcul de la derivee de la dilatation thermique /temps a t=t+TAU/2
C------------------------------------------------------------------
           CALL DERTRA(NYALFA,YALFA,TI12,ALFA1,ALFAV1,TO,TO)
           CALL ZDANUL(EPSTHD,NSTRS1)
           CTH=(ALFAV1*TPOINT*(TI12-TREF))+(ALFA1*TPOINT)
           DO I=1,3
             EPSTHD(I)=CTH
           ENDDO
C
c           DO I=1,3
c             WRITE(6,*) 'sig1(',I,')= ',sig1(I)
c           ENDDO
           CALL XXCREE(TAU,SIG1,EPSV1,VAR1,SIGP2,EVP2,VARP2,EPSTHD,
     &                 DSPT,XMAT,XMAT1,XMAT2,YKX,NKX,NSTRS1,NVARI,
     &                 NCOMAT,NYKX,NNKX,DD,DDV,DDINV,YOG,NYOG,YNU,NYNU,
     &                 MFR1,XCARB,ICARA,IFOURB,2,TI12,TPOINT,TINF,TSUP,
     &                 ITEST,ITHHER,TRUC,NCOURB)
c       Do i=1,3
c         write(6,*) 'sigp2(',i,')=',sigp2(i)
c       enddo
C
           do i=1,nstrs
             sigp2(i) = 0.5d0* ( sigp2(i)+sigp1(i))
             evp2(i) = 0.5d0* ( evp2(i)+evp1(i))
           enddo
           do i=1,nvari
             varp2(i)= 0.5D0 * ( varp2(i)+varp1(i))
           enddo
           t=tau2
           CALL AVANXX(TAU2,SIG,EPSV,VAR,SIG12,EPSV12,VAR12,SIGP2,EVP2,
     &                 VARP2,NSTRS1,NVARI)
c       Do i=1,6
c         write(6,*) 'sig12(',i,')=',sig12(i)
c       enddo
c       Do i=1,6
c         write(6,*) 'epsv12(',i,')=',epsv12(i)
c       enddo
c       Do i=1,7
c         write(6,*) 'var12(',i,')=',var12(i)
c       enddo
           if (tau2.ne.t) then
             tau=2.d0*tau2
             goto 143
           endif
*       if (ib.eq.1.and.igau.eq.1)
*     $ write(6,*)'SIg1(1) SIg2 sig3',sig12(1),sig12(2),sig12(3)
           aap = MAX(ABS(SIG12(1)-SIG(1)),ABS(SIG12(2)-SIG(2)))
           aap = max ( aap,ABS (Sig12(3)-SIG(3)))
           IF ( aap . gt . XMAX * 5.) THEN
*             write(6,*)'SIg12(1) SIg12 SI1g3',sig12(1),sig12(2),sig12(3)
             rap = aap / xmax
             TAU= TAU / rap
*             if (ib.eq.1.and.igau.eq.1)write(6,*)'rap tau'
             go  to 143
*             do I=1,nstrs
*               sig1(i)=XMAX*100.
*               sigf(I)=xmax*200.
*             enddo
*             go to 250
           endif
C
           CALL XXCREE(TAU,SIG12,EPSV12,VAR12,SIGP3,EVP3,VARP3,EPSTHD,
     &                 DSPT,XMAT,XMAT1,XMAT2,YKX,NKX,NSTRS1,NVARI,
     &                 NCOMAT,NYKX,NNKX,DD,DDV,DDINV,YOG,NYOG,YNU,NYNU,
     &                 MFR1,XCARB,ICARA,IFOURB,2,TI12,TPOINT,TINF,TSUP,
     &                 ITEST,ITHHER,TRUC,NCOURB)
           CALL AVANXX(TAU2,SIG12,EPSV12,VAR12,SIG13,EPSV13,VAR13,
     &                 SIGP3,EVP3,VARP3,NSTRS1,NVARI)
           DELTAT=TPOINT*TAU
           TI1=TI0+DELTAT
C Initialisation du tableau XMAT(NCOMAT) a l'instant t=t+TAU
C ou le materiau est a la temperature TI1 comprise dans [TINF,TSUP]
C-----------------------------------------------------------
           CALL INITT(YOG,NYOG,YNU,NYNU,YN,NYN,YM,NYM,YKK,NYKK,
     &                YALFA1,NYALF1,YBETA1,NYBET1,YR,NYR,YA,NYA,YALFA,
     &                NYALFA,YRHO,NYRHO,SIGY,NSIGY,YKX,NYKX,NKX,NNKX,
     &                XMAT,NCOMAT,TI1,TINF,TSUP,ITEST,TRUC,NCOURB)
           IF (ITHHER.EQ.2) THEN
C********** materiau dependant de la temperature **********************
C---------- Initialisation du tableau XMAT1(NCOMAT) a T=TINF
             CALL INITT(YOG,NYOG,YNU,NYNU,YN,NYN,YM,NYM,YKK,NYKK,
     &                  YALFA1,NYALF1,YBETA1,NYBET1,YR,NYR,YA,NYA,
     &                  YALFA,NYALFA,YRHO,NYRHO,SIGY,NSIGY,YKX,NYKX,
     &                  NKX,NNKX,XMAT1,NCOMAT,TINF,TO,TO,ITEST,TRUC,
     &                  NCOURB)
C---------- Initialisation du tableau XMAT2(NCOMAT) a T=TSUP
             CALL INITT(YOG,NYOG,YNU,NYNU,YN,NYN,YM,NYM,YKK,NYKK,
     &                  YALFA1,NYALF1,YBETA1,NYBET1,YR,NYR,YA,NYA,
     &                  YALFA,NYALFA,YRHO,NYRHO,SIGY,NSIGY,YKX,NYKX,
     &                  NKX,NNKX,XMAT2,NCOMAT,TSUP,TO,TO,ITEST,TRUC,
     &                  NCOURB)
C**********************************************************************
           ENDIF
C------------------------------------------------------------------
C Calcul de la derivee de la dilatation thermique /temps a t=t+TAU
C------------------------------------------------------------------
           CALL DERTRA(NYALFA,YALFA,TI1,ALFA1,ALFAV1,TO,TO)
           CALL ZDANUL(EPSTHD,NSTRS1)
           CTH=(ALFAV1*TPOINT*(TI1-TREF))+(ALFA1*TPOINT)
           DO I=1,3
             EPSTHD(I)=CTH
           ENDDO
C
           CALL XXCREE(TAU,SIG13,EPSV13,VAR13,SIGP4,EVP4,VARP4,EPSTHD,
     &                 DSPT,XMAT,XMAT1,XMAT2,YKX,NKX,NSTRS1,NVARI,
     &                 NCOMAT,NYKX,NNKX,DD,DDV,DDINV,YOG,NYOG,YNU,NYNU,
     &                 MFR1,XCARB,ICARA,IFOURB,2,TI1,TPOINT,TINF,TSUP,
     &                 ITEST,ITHHER,TRUC,NCOURB)
           do i=1,nstrs
             sigp4(i)= 0.5D0*( sigp4(i)+sigp3(i))
             evp4(i)  =0.5D0*(  evp4(i)+ evp3(i))
           enddo
           do i=1,nvari
             varp4(i)=0.5D0 * ( varp4(i)+varp3(i))
           enddo
           call avanxx(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,SIGP4,
     &                 EVP4,VARP4,NSTRS1,NVARI)
           CALL XXCREE(TAU,SIG13,EPSV13,VAR13,SIGP4,EVP4,VARP4,EPSTHD,
     &                 DSPT,XMAT,XMAT1,XMAT2,YKX,NKX,NSTRS1,NVARI,
     &                 NCOMAT,NYKX,NNKX,DD,DDV,DDINV,YOG,NYOG,YNU,NYNU,
     &                 MFR1,XCARB,ICARA,IFOURB,2,TI1,TPOINT,TINF,TSUP,
     &                 ITEST,ITHHER,TRUC,NCOURB)
           DO I=1,NSTRS1
             EVP2(I) = (EVP1(I)+EVP4(I))/6.D0+EVP3(I)*2.D0/3.D0
             SIGP2(I)=(SIGP1(I)+SIGP4(I))/6.D0+SIGP3(I)*2.D0/3.D0
           ENDDO
           DO I=1,NVARI
             VARP2(I)=(VARP1(I)+VARP4(I))/6.D0+VARP3(I)*2.D0/3.D0
           ENDDO
           T=TAU
           CALL AVANXX(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,SIGP2,EVP2,
     &                 VARP2,NSTRS1,NVARI)
c       Do i=1,6
c         write(6,*) 'sig1(',i,')=',sig1(i)
c       enddo
c       Do i=1,6
c         write(6,*) 'epsv1(',i,')=',epsv1(i)
c       enddo
c       Do i=1,3
c         write(6,*) 'var1(',i,')=',var1(i)
c       enddo
c           IF (TAU.NE.T) GOTO 143
C---------
         ELSE IF (INPLAS.EQ.142) THEN
C
 144       TAU2=0.5D0*TAU
C           write(6,*) 'tau=',tau
C           write(6,*) 'NSSINC=',NSSINC
*      if (ib.eq.1.and.igau.eq.1)write (6,*) ' tau '  ,tau
           CALL AVANX1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,SIGP1,EVP1,
     &                 VARP1,NSTRS1,NVARI)
c       Do i=2,4
c         write(6,*) 'varp1(',i,')=',varp1(i)
c       enddo
c       Do i=1,6
c         write(6,*) 'sig1(',i,')=',sig1(i)
c       enddo
c       Do i=1,6
c         write(6,*) 'epsv1(',i,')=',epsv1(i)
c       enddo
c       Do i=1,7
c         write(6,*) 'var1(',i,')=',var1(i)
c       enddo
           aap = MAX(ABS(SIG1(1)-SIG(1)),ABS(SIG1(2)-SIG(2)))
           aap = max ( aap,ABS (Sig1(3)-SIG(3)))
           IF ( aap . gt . XMAX * 5.) THEN
*       write(6,*)'sig1(1) sig2 sig3',sig1(1),sig1(2),sig1(3)
             rap = aap / xmax
             TAU= TAU / rap
             go  to 144
*             do I=1,nstrs
*             sig1(i)=XMAX*100.
*             sigf(I)=xmax*200.
*             enddo
*             go to 250
           endif
C
C-----------------------------------------------------------
C Mise a jour eventuelle de la nouvelle temperature TI1
C-----------------------------------------------------------
C
           DELTAT=TPOINT*TAU2
           TI12=TI0+DELTAT
C Initialisation du tableau XMAT(NCOMAT) a l'instant t=t+TAU/2
C ou le materiau est a la temperature TI12 comprise dans [TINF,TSUP]
C-----------------------------------------------------------
           CALL INITT1(YOG1,NYOG1,YNU1,NYNU1,YN1,NYN1,YM1,NYM1,YKK1,
     &                NYKK1,YALF2,NYALF2,YBET2,NYBET2,YR1,NYR1,YA1,NYA1,
     &                YQ1,NYQ1,YALFT1,NYALFT1,YRHO1,NYRHO1,SIGY1,NSIGY1,
     &                XMAT,NCOMAT,TI12,TINF,TSUP,ITEST)
C
           IF (ITHHER.EQ.2) THEN
C********** materiau dependant de la temperature **********************
C---------- Initialisation du tableauXMAT1F(NCOMAT) a T=TINF
             CALL INITT1(YOG1,NYOG1,YNU1,NYNU1,YN1,NYN1,YM1,NYM1,YKK1,
     &                  NYKK1,YALF2,NYALF2,YBET2,NYBET2,YR1,NYR1,YA1,
     &                  NYA1,YQ1,NYQ1,YALFT1,NYALFT1,YRHO1,NYRHO1,
     &                  SIGY1,NSIGY1,XMAT1,NCOMAT,TINF,TO,TO,ITEST)
C---------- Initialisation du tableauXMAT2P(NCOMAT) a T=TSUP
             CALL INITT1(YOG1,NYOG1,YNU1,NYNU1,YN1,NYN1,YM1,NYM1,YKK1,
     &                  NYKK1,YALF2,NYALF2,YBET2,NYBET2,YR1,NYR1,YA1,
     &                  NYA1,YQ1,NYQ1,YALFT1,NYALFT1,YRHO1,NYRHO1,
     &                  SIGY1,NSIGY1,XMAT2,NCOMAT,TSUP,TO,TO,ITEST)
C**********************************************************************
           ENDIF
C------------------------------------------------------------------
C Calcul de la derivee de la dilatation thermique /temps a t=t+TAU/2
C------------------------------------------------------------------
           CALL DERTRA(NYALFT1,YALFT1,TI12,ALFA1,ALFAV1,TO,TO)
           CALL ZDANUL(EPSTHD,NSTRS1)
           CTH=(ALFAV1*TPOINT*(TI12-TREF))+(ALFA1*TPOINT)
           DO I=1,3
             EPSTHD(I)=CTH
           ENDDO
C
c           DO I=1,4
c             WRITE(6,*) 'VAR1(',I,')= ',VAR1(I)
c           ENDDO
c           DO I=1,6
c             WRITE(6,*) 'SIG1(',I,')= ',SIG1(I)
c           ENDDO
           CALL XXCRE1(TAU,SIG1,EPSV1,VAR1,SIGP2,EVP2,VARP2,EPSTHD,
     &                 DSPT,XMAT,XMAT1,XMAT2,NSTRS1,NVARI,NCOMAT,DD,
     &                 DDV,DDINV,DDINVp,YOG1,NYOG1,YNU1,NYNU1,MFR1,
     &                 XCARB,ICARA,IFOURB,2,TI12,TPOINT,TINF,TSUP,
     &                 ITEST,ITHHER,VART1,IB,IGAU,kerre)
c           DO I=1,4
c             WRITE(6,*) 'VAR1 ap(',I,')= ',VAR1(I)
c           ENDDO
c           DO I=1,4
c             WRITE(6,*) 'VARP2(',I,')= ',VARP2(I)
c           ENDDO
c           DO I=1,6
c             WRITE(6,*) 'SIGP2(',I,')= ',SIGP2(I)
c           ENDDO
c           DO I=1,6
c             WRITE(6,*) 'EVP2(',I,')= ',EVP2(I)
c           ENDDO
C
           do i=1,nstrs
             sigp2(i) = 0.5d0* ( sigp2(i)+sigp1(i))
             evp2(i) = 0.5d0* ( evp2(i)+evp1(i))
           enddo
           do i=1,nvari
             varp2(i)= 0.5D0 * ( varp2(i)+varp1(i))
           enddo
           t=tau2
           CALL AVANX1(TAU2,SIG,EPSV,VAR,SIG12,EPSV12,VAR12,SIGP2,EVP2,
     &                 VARP2,NSTRS1,NVARI)
c       Do i=1,6
c         write(6,*) 'sig12(',i,')=',sig12(i)
c       enddo
c       Do i=1,6
c         write(6,*) 'epsv12(',i,')=',epsv12(i)
c       enddo
c       Do i=1,7
c         write(6,*) 'var12(',i,')=',var12(i)
c       enddo
           if (tau2.ne.t) then
             tau=2.d0*tau2
             goto 144
           endif
*       if (ib.eq.1.and.igau.eq.1)
*     $ write(6,*)'SIg1(1) SIg2 sig3',sig12(1),sig12(2),sig12(3)
           aap = MAX(ABS(SIG12(1)-SIG(1)),ABS(SIG12(2)-SIG(2)))
           aap = max ( aap,ABS (Sig12(3)-SIG(3)))
           IF ( aap . gt . XMAX * 5.) THEN
*             write(6,*)'SIg12(1) SIg12 SI1g3',sig12(1),sig12(2),sig12(3)
             rap = aap / xmax
             TAU= TAU / rap
*             if (ib.eq.1.and.igau.eq.1)write(6,*)'rap tau'
             go  to 144
*             do I=1,nstrs
*               sig1(i)=XMAX*100.
*               sigf(I)=xmax*200.
*             enddo
*             go to 250
           endif
C
           CALL XXCRE1(TAU,SIG12,EPSV12,VAR12,SIGP3,EVP3,VARP3,EPSTHD,
     &                 DSPT,XMAT,XMAT1,XMAT2,NSTRS1,NVARI,NCOMAT,DD,
     &                 DDV,DDINV,DDINVp,YOG1,NYOG1,YNU1,NYNU1,MFR1,
     &                 XCARB,ICARA,IFOURB,2,TI12,TPOINT,TINF,TSUP,
     &                 ITEST,ITHHER,VART1,IB,IGAU,kerre)
           CALL AVANX1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPSV13,VAR13,
     &                 SIGP3,EVP3,VARP3,NSTRS1,NVARI)
           DELTAT=TPOINT*TAU
           TI1=TI0+DELTAT
 
C Initialisation du tableau XMAT(NCOMAT) a l'instant t=t+TAU
C ou le materiau est a la temperature TI1 comprise dans [TINF,TSUP]
C-----------------------------------------------------------
           CALL INITT1(YOG1,NYOG1,YNU1,NYNU1,YN1,NYN1,YM1,NYM1,YKK1,
     &                 NYKK1,YALF2,NYALF2,YBET2,NYBET2,YR1,NYR1,YA1,
     &                 NYA1,YQ1,NYQ1,YALFT1,NYALFT1,YRHO1,NYRHO1,SIGY1,
     &                 NSIGY1,XMAT,NCOMAT,TI1,TINF,TSUP,ITEST)
           IF (ITHHER.EQ.2) THEN
C********** materiau dependant de la temperature **********************
C---------- Initialisation du tableau XMAT1(NCOMAT) a T=TINF
             CALL INITT1(YOG1,NYOG1,YNU1,NYNU1,YN1,NYN1,YM1,NYM1,YKK1,
     &                   NYKK1,YALF2,NYALF2,YBET2,NYBET2,YR1,NYR1,YA1,
     &                   NYA1,YQ1,NYQ1,YALFT1,NYALFT1,YRHO1,NYRHO1,
     &                   SIGY1,NSIGY1,XMAT1,NCOMAT,TINF,TO,TO,ITEST)
C---------- Initialisation du tableau XMAT2(NCOMAT) a T=TSUP
             CALL INITT1(YOG1,NYOG1,YNU1,NYNU1,YN1,NYN1,YM1,NYM1,YKK1,
     &                   NYKK1,YALF2,NYALF2,YBET2,NYBET2,YR1,NYR1,YA1,
     &                   NYA1,YQ1,NYQ1,YALFT1,NYALFT1,YRHO1,NYRHO1,
     &                   SIGY1,NSIGY1,XMAT2,NCOMAT,TSUP,TO,TO,ITEST)
C**********************************************************************
           ENDIF
C------------------------------------------------------------------
C Calcul de la derivee de la dilatation thermique /temps a t=t+TAU
C------------------------------------------------------------------
           CALL DERTRA(NYALFT1,YALFT1,TI1,ALFA1,ALFAV1,TO,TO)
           CALL ZDANUL(EPSTHD,NSTRS1)
           CTH=(ALFAV1*TPOINT*(TI1-TREF))+(ALFA1*TPOINT)
           DO I=1,3
             EPSTHD(I)=CTH
           ENDDO
C
           CALL XXCRE1(TAU,SIG13,EPSV13,VAR13,SIGP4,EVP4,VARP4,EPSTHD,
     &                 DSPT,XMAT,XMAT1,XMAT2,NSTRS1,NVARI,NCOMAT,DD,
     &                 DDV,DDINV,DDINVp,YOG1,NYOG1,YNU1,NYNU1,MFR1,
     &                 XCARB,ICARA,IFOURB,2,TI1,TPOINT,TINF,TSUP,
     &                 ITEST,ITHHER,VART1,IB,IGAU,kerre)
           do i=1,nstrs
             sigp4(i)= 0.5D0*( sigp4(i)+sigp3(i))
             evp4(i)  =0.5D0*(  evp4(i)+ evp3(i))
           enddo
           do i=1,nvari
             varp4(i)=0.5D0 * ( varp4(i)+varp3(i))
           enddo
           call AVANX1(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,SIGP4,
     &                 EVP4,VARP4,NSTRS1,NVARI)
           CALL XXCRE1(TAU,SIG13,EPSV13,VAR13,SIGP4,EVP4,VARP4,EPSTHD,
     &                 DSPT,XMAT,XMAT1,XMAT2,NSTRS1,NVARI,NCOMAT,DD,
     &                 DDV,DDINV,DDINVp,YOG1,NYOG1,YNU1,NYNU1,MFR1,
     &                 XCARB,ICARA,IFOURB,2,TI1,TPOINT,TINF,TSUP,
     &                 ITEST,ITHHER,VART,IB,IGAU,kerre)
           DO I=1,NSTRS1
             EVP2(I) = (EVP1(I)+EVP4(I))/6.D0+EVP3(I)*2.D0/3.D0
             SIGP2(I)=(SIGP1(I)+SIGP4(I))/6.D0+SIGP3(I)*2.D0/3.D0
           ENDDO
           DO I=1,NVARI
             VARP2(I)=(VARP1(I)+VARP4(I))/6.D0+VARP3(I)*2.D0/3.D0
           ENDDO
           T=TAU
           CALL AVANX1(TAU,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,SIGP2,EVP2,
     &                 VARP2,NSTRS1,NVARI)
c       Do i=1,6
c         write(6,*) 'sig1(',i,')=',sig1(i)
c       enddo
c       Do i=1,6
c         write(6,*) 'epsv1(',i,')=',epsv1(i)
c       enddo
c       Do i=1,7
c         write(6,*) 'var1(',i,')=',var1(i)
c       enddo
c           IF (TAU.NE.T) GOTO 144
C---------
         ELSE IF (INPLAS.EQ.165) THEN
           tau2=tau*0.5d0
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,ZMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           DO I=1,NCOMAT
             ZMAT(I)  = 0.5D0*(XMAT(I)+XMAT0(I))
           ENDDO
           CALL INCRE9(SIG1,VAR1,ZMAT,NSTRS1,EVP2,VARP2)
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
           ENDDO
           DO I=1,NVARI
            VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
           ENDDO
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
     &                 VARP2,ZMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE9(SIG12,VAR12,ZMAT,NSTRS1,EVP3,VARP3)
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,ZMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           CALL INCRE9(SIG13,VAR13,ZMAT,NSTRS1,EVP4,VARP4)
           DO I=1,NSTRS1
             EVP4(I)  = 0.5d0*(EVP3(I)+EVP4(I))
           ENDDO
           DO I=1,NVARI
             VARP4(I) = 0.5d0*(VARP3(I)+VARP4(I))
           ENDDO
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARF,DSPT,
     &                 EVP4,VARP4,ZMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
           DO I=1,NCOMAT
             ZMAT(I)  = XMAT(I)
           ENDDO
           CALL INCRE9(SIGF,VARF,ZMAT,NSTRS1,EVP4,VARP4)
           DO I=1,NSTRS1
             EVP2(I)  = (EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           ENDDO
           DO I=1,NVARI
             VARP2(I) = (VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
           ENDDO
           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                 VARP2,ZMAT,NSTRS1,NVARI,IFOURB,INPLAS,MFR1)
C---------
         ENDIF
         GOTO 250
C _____________________________________________________________________
150    CONTINUE
C----------------------------------------------------------------------
C  CALCULATIONS FOR GENERALISED STRESS/STRAIN FORMULATIONS
C----------------------------------------------------------------------
         IF ((INPLAS.GE.19.AND.INPLAS.LE.24).OR.INPLAS.EQ.61) THEN
           tau2=tau*0.5d0
           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
           CALL INCRE3(TAU2,SIG1,EPSV1,VAR1,XMAT,EVP2,VARP2,ALFA,
     &                 NSTRS1,NVARI,INPLAS,NCOMAT)
           DO I=1,NSTRS1,1
             EVP2(I) = 0.5d0*(EVP1(I)+EVP2(I))
           ENDDO
           IF (INPLAS.EQ.24) THEN
             DO I=1,NVARI
               VARP2(I)=0.5D0*(VARP1(I)+VARP2(I))
             ENDDO
           ELSE
             DO I=1,2*NSTRS1+2
               VARP2(I) = 0.5d0*(VARP1(I)+VARP2(I))
             ENDDO
             DO I=2*NSTRS1+4,NVARI
               VARP2(I) = 0.5d0*(VARP1(I)+VARP2(I))
             ENDDO
           ENDIF
           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
           CALL INCRE3(TAU2,SIG12,EPSV12,VAR12,XMAT,EVP3,VARP3,ALFA,
     &                 NSTRS1,NVARI,INPLAS,NCOMAT)
           CALL ADVAN2(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
     &                 NCOMAT)
           CALL INCRE3(TAU2,SIG13,EPSV13,VAR13,XMAT,EVP4,VARP4,ALFA,
     &                 NSTRS1,NVARI,INPLAS,NCOMAT)
           DO I=1,NSTRS1
              EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
           enddo
           IF (INPLAS.EQ.24) THEN
             DO I=1,NVARI
               VARP4(I)=0.5D0*(VARP3(I)+VARP4(I))
             enddo
           ELSE
             DO I=1,2*NSTRS1+2
              VARP4(I) = 0.5d0*(VARP3(I)+VARP4(I))
             enddo
             DO I=2*NSTRS1+4,NVARI
               VARP4(I) = 0.5d0*(VARP3(I)+VARP4(I))
             enddo
           ENDIF
           CALL ADVAN2(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
     &                 NCOMAT)
           CALL INCRE3(TAU2,SIGf,EPinf,VARf,XMAT,EVP4,VARP4,ALFA,
     &                 NSTRS1,NVARI,INPLAS,NCOMAT)
           DO I=1,NSTRS1,1
             EVP2(I) =(EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           enddo
           IF (INPLAS.EQ.24) THEN
             DO I=1,NVARI
               VARP2(I)=(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
             enddo
           ELSE
             DO I=1,2*NSTRS1+2
               VARP2(I) =(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
             enddo
             DO I=2*NSTRS1+4,NVARI
               VARP2(I) =(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
             enddo
           endif
           CALL ADVAN2(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                 VARp2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
C---------
         ELSE IF (INPLAS.EQ.25) THEN
           tau2=tau*0.5d0
           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
           CALL INCRE4(SIG1,VAR1,EVP2,VARP2,XMAT,ALFA,NSTRS1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP2(I) = 0.5d0*(EVP1(I)+EVP2(I))
           ENDDO
           DO I=1,NVARI
             VARP2(I) = 0.5d0*(VARP1(I)+VARP2(I))
           ENDDO
           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
           CALL INCRE4(SIG12,VAR12,EVP3,VARP3,XMAT,ALFA,NSTRS1,NVARI,
     &                 NCOMAT)
           CALL ADVAN2(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
     &                 NCOMAT)
           CALL INCRE4(SIG13,VAR13,EVP4,VARP4,XMAT,ALFA,NSTRS1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
           enddo
           DO I=1,NVARI
             VARP4(I) = 0.5d0*(VARP3(I)+VARP4(I))
           enddo
           CALL ADVAN2(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
     &                 NCOMAT)
           CALL INCRE4(SIGf,VARf,EVP4,VARP4,XMAT,ALFA,NSTRS1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP2(I) =(EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           enddo
           DO I=1,NVARI
             VARP2(I) =(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
           enddo
           CALL ADVAN2(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
C---------
         ELSE IF (INPLAS.EQ.76) THEN
           CALL ADVAN2(TAU,SIG,EPSV,VAR,SIGf,EPinf,VARf,DSPT,
     &                 EVP1,VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
     &                 NCOMAT)
           CALL INCRA4(SIGf,VARf,EVP2,VARP2,XMAT,ALFA,NSTRS1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
           ENDDO
           DO I=1,NVARI
             VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
           ENDDO
           CALL ADVAN2(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
*           tau2=tau*0.5d0
*           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,
*     &                 EVP1,VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
*     &                 NCOMAT)
*           CALL INCRA4(SIG1,VAR1,EVP2,VARP2,XMAT,ALFA,NSTRS1,NVARI,
*     &                 NCOMAT)
*           DO I=1,NSTRS1
*             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
*           ENDDO
*           DO I=1,NVARI
*             VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
*           ENDDO
*           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
*     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
*           CALL INCRA4(SIG12,VAR12,EVP3,VARP3,XMAT,ALFA,NSTRS1,NVARI,
*     &                 NCOMAT)
*           CALL ADVAN2(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
*     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
*     &                 NCOMAT)
*           CALL INCRA4(SIG13,VAR13,EVP4,VARP4,XMAT,ALFA,NSTRS1,NVARI,
*     &                 NCOMAT)
*           DO I=1,NSTRS1
*              EVP4(I) = 0.5D0*(EVP3(I)+EVP4(I))
*           enddo
*           DO I=1,NVARI
*             VARP4(I) = 0.5D0*(VARP3(I)+VARP4(I))
*           enddo
*           CALL ADVAN2(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
*     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
*     &                 NCOMAT)
*           DO I=1,NSTRS1
*             EVP2(I) = (EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
*           enddo
*           DO I=1,NVARI
*             VARP2(I) = (VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
*           enddo
*           CALL ADVAN2(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
*     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
C---------
         ELSE IF (INPLAS.EQ.77) THEN
           CALL ADVAN2(TAU,SIG,EPSV,VAR,SIGf,EPinf,VARf,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,ncomat)
           CALL INCRB4(SIGf,VARf,EVP2,VARP2,XMAT,ALFA,NSTRS1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP2(I) = 0.5d0*(EVP1(I)+EVP2(I))
           enddo
           DO I=1,NVARI
             VARP2(I) = 0.5d0*(VARP1(I)+VARP2(I))
           enddo
           CALL ADVAN2(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
*           tau2=tau*0.5d0
*           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
*     &                 VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
*           CALL INCRB4(SIG1,VAR1,EVP2,VARP2,XMAT,ALFA,NSTRS1,NVARI,
*     &                 NCOMAT)
*           DO I=1,NSTRS1
*             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
*           ENDDO
*           DO I=1,NVARI
*             VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
*           ENDDO
*           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
*     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
*           CALL INCRB4(SIG12,VAR12,EVP3,VARP3,XMAT,ALFA,NSTRS1,NVARI,
*     &                 NCOMAT)
*           CALL ADVAN2(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
*     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
*     &                 NCOMAT)
*           CALL INCRB4(SIG13,VAR13,EVP4,VARP4,XMAT,ALFA,NSTRS1,NVARI,
*     &                 NCOMAT)
*           DO I=1,NSTRS1
*             EVP4(I) = 0.5D0*(EVP3(I)+EVP4(I))
*           enddo
*           DO I=1,NVARI
*             VARP4(I) = 0.5D0*(VARP3(I)+VARP4(I))
*           enddo
*           CALL ADVAN2(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
*     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
*     &                 NCOMAT)
*           CALL INCRB4(SIGf,VARf,EVP4,VARP4,XMAT,ALFA,NSTRS1,NVARI,
*     &                 NCOMAT)
*           DO I=1,NSTRS1
*             EVP2(I) =(EVP3(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
*           enddo
*           DO I=1,NVARI
*             VARP2(I) =(VARP3(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
*           enddo
*           CALL ADVAN2(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
*     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
C---------
         ELSE IF (INPLAS.EQ.53) THEN
           tau2=tau*0.5d0
           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,
     &                 EVP1,VARP1,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
     &                 NCOMAT)
           CALL INCRE6(SIG1,VAR1,EVP2,VARP2,XMAT,ALFA,NSTRS1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
           ENDDO
           DO I=1,NVARI
             VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
           ENDDO
           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
           CALL INCRE6(SIG12,VAR12,EVP3,VARP3,XMAT,ALFA,NSTRS1,NVARI,
     &                 NCOMAT)
           CALL ADVAN2(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
     &                 NCOMAT)
           CALL INCRE6(SIG13,VAR13,EVP4,VARP4,XMAT,ALFA,NSTRS1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP4(I) = 0.5D0*(EVP3(I)+EVP4(I))
           enddo
           DO I=1,NVARI
             VARP4(I) = 0.5D0*(VARP3(I)+VARP4(I))
           enddo
           CALL ADVAN2(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,
     &                 EVP4,VARP4,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,
     &                 NCOMAT)
           CALL INCRE6(SIGf,VARf,EVP4,VARP4,XMAT,ALFA,NSTRS1,NVARI,
     &                 NCOMAT)
           DO I=1,NSTRS1
             EVP2(I) =(EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           enddo
           DO I=1,NVARI
             VARP2(I) =(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
           enddo
           CALL ADVAN2(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,
     &                 VARP2,XMAT,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
C---------
         ENDIF
         GOTO 250
C
C CAS D'UN TUYAU FISSURE EN FLUAGE SUIVANT LA LOI "NORTON"
C
 210     CONTINUE
         IF (INPLAS.NE.19) GOTO 999
         tau2=tau*0.5d0
         CALL TUFADV(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,VARP1,
     &               XMAT,XCARB,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
         CALL TUFINC(TAU2,SIG1,EPSV1,VAR1,XMAT,XCARB,EVP2,VARP2,NSTRS1,
     &               NVARI,INPLAS,NCOMAT,KERREU1)
         DO I=1,6
           EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
         ENDDO
         DO I=1,4
           VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
         ENDDO
         DO I=6,NVARI
           VARP2(I) = 0.5D0*(VARP1(I)+VARP2(I))
         ENDDO
         CALL TUFADV(TAU2,SIG,EPSV,VAR,SIG12,EPsv12,VAR12,DSPT,EVP2,
     &               VARP2,XMAT,XCARB,NSTRS1,NVARI,IFOURB,INPLAS,
     &               NCOMAT)
         CALL TUFINC(TAU2,SIG12,EPSV12,VAR12,XMAT,XCARB,EVP3,VARP3,
     &               NSTRS1,NVARI,INPLAS,NCOMAT,KERREU1)
         CALL TUFADV(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &               EVP3,VARP3,XMAT,XCARB,NSTRS1,NVARI,IFOURB,INPLAS,
     &               NCOMAT)
         CALL TUFINC(TAU2,SIG13,EPSV13,VAR13,XMAT,XCARB,EVP4,VARP4,
     &               NSTRS1,NVARI,INPLAS,NCOMAT,KERREU1)
         DO I=1,6
           EVP4(I) = 0.5D0*(EVP3(I)+EVP4(I))
         enddo
         DO I=1,4
           VARP4(I) = 0.5D0*(VARP3(I)+VARP4(I))
         enddo
         DO I=6,NVARI
           VARP4(I) = 0.5D0*(VARP3(I)+VARP4(I))
         enddo
         CALL TUFADV(TAU2,SIG12,EPSV12,VAR12,SIGf,EPinf,VARf,DSPT,EVP4,
     &               VARP4,XMAT,XCARB,NSTRS1,NVARI,IFOURB,INPLAS,
     &               NCOMAT)
         CALL TUFINC(TAU2,SIGf,EPinf,VARf,XMAT,XCARB,EVP4,VARP4,NSTRS1,
     &               NVARI,INPLAS,NCOMAT,KERREU1)
         DO I=1,6
           EVP2(I) =(EVP3(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
         enddo
         DO I=1,4
           VARP2(I) =(VARP3(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
         enddo
         DO I=6,NVARI
           VARP2(I) =(VARP3(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
         enddo
         CALL TUFADV(TAU,SIG,EPSV,VAR,SIG1,EPsv1,VAR1,DSPT,EVP2,VARP2,
     &               XMAT,XCARB,NSTRS1,NVARI,IFOURB,INPLAS,NCOMAT)
C
C----------------------------------------------------------------------
 250   CONTINUE
C ---------------------------------------------------------------------
C   CALCUL DU RAPPORT : ERREUR CALCULEE / ERREUR ADMISE
C ---------------------------------------------------------------------
****       CALL ZERO(XX  ,8    ,1)      inutile?
       DO 2501 I=1,NSTRS1     
         CALL VERREE(SIG1(I),ILOG)
         IF (.NOT.ILOG) GOTO 2502
 
C        Copie de SIGF(I) dans XTEST sinon ARGUMENT ELEMENT DE SEGMENT
         CALL VERREE(SIGF(I),ILOG)
         IF (.NOT.ILOG) GOTO 259
         XX(I) = SIGF(I)-SIG1(I)
 2501  CONTINUE
 
 2502  CONTINUE
       ITAIL= MIN(I,NSTRS1)
       RA   = SQRT(PROCON(XX,XX,ITAIL))/(ERRABS)
c
       IF (MFR1.EQ.17 .AND. INPLAS.EQ.19) THEN
         RA=SQRT(XX(1)**2 + XX(6)**2)/ERRABS
       ENDIF
       SQRA = SQRT(RA)
C
       IF (INPLAS.EQ.29) THEN
         RD = 0.D0
         IF (VARF(3).NE.0.D0.AND.VAR1(3).GT.1.0D-5) THEN
           RD = (VARF(3)-VAR1(3))/VARF(3)
           RD = ABS(RD)/1.D-4
         ENDIF
         RA =  MAX(RA,RD)
         SQRA = SQRT(RA)
       ELSEIF (INPLAS.EQ.142) THEN
         RD = 0.D0
         IF (VARF(8).NE.0.D0.AND.VAR1(8).GT.1.0D-5) THEN
           RD = (VARF(8)-VAR1(8))/VARF(8)
           RD = ABS(RD)/1.D-4
         ENDIF
         RA =  MAX(RA,RD)
         SQRA = SQRT(RA)
       ELSEIF (INPLAS.EQ.25) THEN
C        Convergence variables internes ONERA
         RD = 0.D0
C        Uniquement si multiplicateur plastique croit sensiblement
C        Cas des ecrouissages cinematiques, homogenes a une contrainte
C        ERRABS : precision sur la contrainte
C        PRELOC : precision relative
C        XN     : ordre de grandeur de la contrainte
         XN = ERRABS/PRELOC
         DO 251 I=1,2*NSTRS1
           XI = MAX(VARF(I),XN)
           VR = ABS(VARF(I)-VAR1(I))/XI
           IF (VR.GT.RD) RD = VR
 251     CONTINUE
C
C        Cas de l'ecrouissage isotrope R
         I = 4*NSTRS1+2
         XI = MAX(VARF(I),XN)
         VR = ABS(VARF(I)-VAR1(I))/XI
         IF (VR.GT.RD) RD = VR
C
C        Cas de la variable QQQ
         I = 4*NSTRS1+3
         XI = MAX(VARF(I),XN)
         VR = ABS(VARF(I)-VAR1(I))/XI
         IF (VR.GT.RD) RD=VR
C
         RD = RD/ERRABS
C
C        Autres variables internes : homognes a des deformations
C        ERRABS2 : precision sur les deformations (1.D11~module d'Young)
C        XN2     : ordre de grandeur des deformations
         RD2 = 0.D0
         ERRABS2 = MAX(ERRABS/1.D11,1.D-12)
         XN2 = ERRABS2/PRELOC
C
         DO 252 I=2*NSTRS1+1,4*NSTRS1
           XI2 = MAX(VARF(I),XN2)
           VR = ABS(VARF(I)-VAR1(I))/XI2
           IF (VR.GT.RD2) RD2 = VR
 252     CONTINUE
C
C        Cas de la deformation plastique cumulee
         I = 4*NSTRS1+1
         XI2 = MAX(VARF(I),XN2)
         VR = ABS(VARF(I)-VAR1(I))/XI2
         IF (VR.GT.RD2) RD2=VR
C
C        Cas de la variable QQ
         I = 4*NSTRS1+4
         XI2 = MAX(VARF(I),XN2)
         VR = ABS(VARF(I)-VAR1(I))/XI2
         IF (VR.GT.RD2) RD2=VR
C
         RD2 = RD2/ERRABS2
         IF (RD2.GT.RD) RD = RD2
C
         RA =  MAX(RA,RD)
         SQRA = SQRT(RA)
       ELSEIF (INPLAS.EQ.165) THEN
C        Convergence variables internes Chaboche
         RD = 0.D0
C        Les variables internes sont homognes a des deformations
C        ERRABS  : precision sur les contraintes
C        ERRABS2 : precision sur les deformations (1.D11~module d'Young)
C        PRELOC  : precision relative
C        XN2     : ordre de grandeur des deformations
         ERRABS2 = MAX(ERRABS/1.D11,1.D-12)
         XN2 = ERRABS2/PRELOC
         DO 260 I=1,NVARI
             XI2 = MAX(VARF(I),XN2)
             VR = ABS(VARF(I)-VAR1(I))/XI2
             IF (VR.GT.RD) RD = VR
 260     CONTINUE
         RD = RD/ERRABS2
         RA =  MAX(RA,RD)
         SQRA = SQRT(RA)
       ENDIF
C ---------------------------------------------------------------------
C  TEST DE FIN D'ITERATIONS /  MISE A JOUR DE TAU  /OPTION JECHER
C  DIV =7   BORNE = 2
C  SI   SQRA>7   TAU = TAU/7    ET NOUVEL ESSAI
C  SI 2&lt;RA&lt;7*7   ON VISE RA = 1  ET NOUVEL ESSAI
C ------------------------------------------------------------------
       IF (.not.iforce.and.dtlibr ) Then
c petite ruse pour dejouer l'optimisation
                 ra1=ra*1.d0
*       write(6,*) ' ra, div tau taux' , ra, div,tau,taux
         IF ((RA.GT.DIV*DIV).OR.(RA.NE.RA1)) THEN
           TAU = TAU/div
           IF ((INPLAS.EQ.29).OR.(INPLAS.EQ.142)) TAU =MIN(TAU,TAUX)
           DELTAT=TPOINT*TAU
           TI1=TI0+DELTAT
           GOTO 80
         ELSEIF ( RA.GT.(BORNE)) THEN
           TAU = TAU/SQRA
           IF ((INPLAS.EQ.29).OR.(INPLAS.EQ.142)) TAU =MIN(TAU,TAUX)
           DELTAT=TPOINT*TAU
           TI1=TI0+DELTAT
           GOTO 80
         ENDIF
       ENDIF
C ---------------------------------------------------------------------
C       ici  ra &lt; borne      cas JECHER  :
C ---------------------------------------------------------------------
C je n'ai pas trouvé comment jecher = 1 pour moi jamais. TC
       IF (JECHER.EQ.1) THEN
*         write(6,*) ' on passe dans jecher = 1'
         DTT = TAU
         NSSINC = NITERA
         IF ((NSSINC.EQ.1).AND.(RA.EQ.0.0)) GOTO 999
         IF (NITERA.GE.8) GOTO 999
         IF (FAC*SQRA.LT.1.0) THEN
           TAU = TAU*FAC
           DELTAT=TPOINT*TAU
           TI1=TI0+DELTAT
           GOTO 80
         ELSEIF ((SQRA.LT.RMIN).OR.(SQRA.GT.RMAX)) THEN
           TAU = TAU/SQRA
           DELTAT=TPOINT*TAU
           TI1=TI0+DELTAT
           GOTO 80
         ENDIF
C ---------------------------------------------------------------------
C         ici  rmin &lt; sqra &lt; rmax    et  nitera &lt; 8
C         pas de mise @ jour des variables
C ---------------------------------------------------------------------
         GOTO 999
       ENDIF
C ----------------------------------------------------------------------
C  FIN D'ITERATIONS  / MISE A JOUR DES VARIABLES
C  ici RA &lt; BORNE
C  fin des boucles sur tau optimal
C  on avance en temps
C  mise @ jour de SIG etc...
C -------------------------------------------------------------------
       INV = INV + IVTEST
       DO I=1,NSTRS1
         SIG(I)  = SIGF(I)
         EPSV(I) = EPINF(I)
       ENDDO
       DO I=1,NVARI
         VAR(I)  = VARF(I)
       ENDDO
*       if(pasbea.lt.ra) pasbea=ra
       if (iforce) iffo=iffo+1
 
*       IF ( nssinc.eq. 1) dtprem = tau
*       IF ( nssinc.eq. 2) dtseco = tau
C
       IF (INPLAS.EQ.29) THEN
C--------------------------------------------------------------
C Estimation du pas de temps apres la mise a jour des variables
C--------------------------------------------------------------
C
         CALL ESTITO(SIG,NSTRS1,VAR,NVARI,YKX,NYKX,NKX,NNKX,XMAT,
     &               NCOMAT,TI1,TD,TRUC,NCOURB)
 
         IF ((VARF(3).GE.0.96).OR.(TD.LT.1.D0)) THEN
           VARF(3)=1.D0
           TLIFE = DT - (DTLEFT - TAU)
           GOTO 999
         ENDIF
C
       ELSEIF (INPLAS.EQ.142) THEN
C--------------------------------------------------------------
C Estimation du pas de temps apres la mise a jour des variables
C--------------------------------------------------------------
C
         CALL ESTIT1(SIG,NSTRS1,VAR,NVARI,XMAT,NCOMAT,TI1,TD)
 
         IF ((VARF(8).GE.0.96).OR.(TD.LT.1.D0)) THEN
           VARF(8)=1.D0
           TLIFE = DT - (DTLEFT - TAU)
           GOTO 999
         ENDIF
       ENDIF
C
C --------------------------------------------------------------------
C  TEST DE FIN SS INCREMENTS /  MISE A JOUR DE TAU
C  si   SQRA&lt;1/3        TAU = TAU*3
C  si   1/3&lt;SQRA&lt;RMIN   on vise RA = 1
C  si   RMIN&lt;SQRA&lt;RMAX  TAU inchang{
C  si   SQRA>RMAX       on vise RA = 1
C    fin des boucles en ss increments  si tau = dtleft
C --------------------------------------------------------------------
C
       IF ( TAU.LT.DTLEFT ) THEN
*         DTDEUX=TAU
         DTLEFT = DTLEFT - TAU
*         IF (dtlibr) then
         IF ( FAC*SQRA.LT.1.D0) THEN
           TAU=TAU*FAC
         ELSEIF ( (SQRA.LT.RMIN).OR.(SQRA.GT.RMAX) ) THEN
           TAU=TAU/SQRA
         ENDIF
*          else
*            TAU = TAU * R
*          endif
         IF (TAU.GT.DTLEFT) then
           TAU = DTLEFT
         endif
         IF ((INPLAS.EQ.29).OR.(INPLAS.EQ.107).OR.(INPLAS.EQ.142)) THEN
C----------------------------------------------------------------------------
C Mise a jour des temperatures
C TI0 temperature au dedut du pas de sous-incrementation avec TINF&lt;TI0&lt;TSUP
C TI1 temperature a la fin du pas de sous-incrementation
C-----------------------------------------------------------------------------
           TI0=TI1
           DELTAT=TPOINT*TAU
           TI1=TI0+DELTAT
         ENDIF
         IF (INPLAS.EQ.107) THEN
C----------------------------------------------------------------------------
C Mise a jour des densites de fissions
C FII0 densite de fissions au dedut du pas de sous-incrementation
C FII1 densite de fissions a la fin du pas de sous-incrementation
C-----------------------------------------------------------------------------
           FII0=FII1
           DELTAF=FPOINT*TAU
           FII1=FII0+DELTAF
         ENDIF
         GOTO 70
       ENDIF
C
       IF (ABS(TAU-DTLEFT).GT.(TAU/1000.)) THEN
         WRITE ( IOIMP,* ) ' PROBLEME   TAU > DTLEFT '
         KERRE = 223
       ENDIF
 
C-----------------------------------------------------------------------
999    CONTINUE
       IF (MFR1.EQ.3) THEN
         DO 1000 I=1,NSTRS1/2
           SIGF(        I) =SIGF(         I)*THICK
           SIGF(NSTRS1/2+I) =SIGF(NSTRS1/2+ I)*THICK*THICK/6.0
*           DSIGT(        I)=DSIGT(        I)*THICK
*           DSIGT(NSTRS1/2+I)=DSIGT(NSTRS1/2+I)*THICK*THICK/6.0
1000     CONTINUE
       ENDIF
C
C===========================================================
C  RETOUR A LA DEFINITION NORMALE DES DEFORMATIONS
C  A SAVOIR: LES DEFORMATIONS DE CISAILLEMENT SONT
C  DEFINIES PAR DES GAMA.
C  ON MULTIPLIE DONC LES TERMES DE CISAILLEMENT PAR 2.
C  CECI NE CONCERNE PAS LE MODELE VISCO-ENDOMMAGEABLE
C  DE LEMAITRE (INPLAS=29).
C
C  SEULES LES FORMULATIONS SUIVANTES SONT ACCEPTEES PAR CONSTI:
C      MFR1=1  (MASSIF)
C      MFR1=5  (COQUES EPAISSES)
C      MFR1=3  (COQUES MINCES)
C      MFR1=17 (TUYAUX FISSURES)
C      MFR1=31 (BBAR)
C      MFR1=33 (POREUX)
C
       IF ((INPLAS.NE.29).AND.(INPLAS.NE.142)) THEN
C
C  Cas de la formulation massive
C  Les termes de cisaillement apparaissent
C  au delà de la troisieme composante
C
         IF (MFR1.EQ.1.OR.MFR1.EQ.31.OR.MFR1.EQ.33) THEN
           DO 14 I=1,NSTRS1
             A=1.D0
             IF (I.GT.3) A=2.D0
             EPIN0(I)=EPIN0(I)*A
             EPINF(I)=EPINF(I)*A
 14        CONTINUE
C
C  Cas des coques épaisses
C  Les termes de cisaillement apparaissent
C  au delà de la deuxieme composante
C
         ELSE IF (MFR1.EQ.5) THEN
           DO 15 I=1,NSTRS1
             A=1.D0
             IF (I.GT.2) A=2.D0
             EPIN0(I)=EPIN0(I)*A
             EPINF(I)=EPINF(I)*A
 15        CONTINUE
C
C  Cas des coques minces
C  Les termes de cisaillement apparaissent
C  pour la troisieme et la sixieme composante
C  uniquement dans les cas de calculs
C  tridimensionnels ou d'analyse de Fourier
C
         ELSE IF (MFR1.EQ.3) THEN
           IF ((IFOURB.EQ.1).OR.(IFOURB.EQ.2)) THEN
             DO 16 I=1,NSTRS1
               A=1.D0
               IF (I.EQ.3) A=2.D0
               IF (I.EQ.6) A=2.D0
               EPIN0(I)=EPIN0(I)*A
               EPINF(I)=EPINF(I)*A
 16          CONTINUE
           ENDIF
C
C  Reste le cas des tuyaux fissurés (MFR1=17)
C
         ENDIF
       ENDIF
C
C===========================================================
C
 998  RETURN
C
C     On detecte des NaN ou des Inf dans les contraintes
 259  CALL ERREUR(460)
      RETURN
C
      END