Sources Esope | Cast3M

cconst
C CCONST    SOURCE    OF166741  25/11/04    21:15:10     12349          
      SUBROUTINE CCONST(wrk52,wrk53,wrk54,WRK7,WRK8,WRK9,
     1     WRK91,NVARI,NSSINC,INV,IFOUR7,T0,TF,FI0,FIF
     4    ,TLIFE,NCOUR7,IB,IGAU,NBPGAU,KERREU1,iecou,xecou)
 
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(NCOUR7) 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         IFOUR7 = -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         TREFA 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
-INC TECOU
 
      SEGMENT WRK7
        REAL*8 F(NCOUR7,2),W(NCOUR7),TRUC(NCOUR7)
      ENDSEGMENT
 
      SEGMENT WRK8
        REAL*8 DD(NSTRS1,NSTRS1),DDV(NSTRS1,NSTRS1),DDINV(NSTRS1,NSTRS1)
        REAL*8 DDINVp(NSTRS1,NSTRS1)
      ENDSEGMENT
 
      SEGMENT WRK9
        REAL*8 YOG(NYOG),YNU(NYNU),YALFA(NYALFA),YSMAX(NYSMAX)
        REAL*8 YN(NYN),YM(NYM),YKK(NYKK),YALFA1(NYALF1)
        REAL*8 YBETA1(NYBET1),YR(NYR),YA(NYA),YKX(NYKX),YRHO(NYRHO)
        REAL*8 SIGY(NSIGY)
        INTEGER NKX(NNKX)
      ENDSEGMENT
 
      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
 
       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),XM(6)
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)
 
       logical dtlibr,iforce,ilog
 
       CALL ZERO(SIG1 ,8    ,1)
 
       NSTRS1 = iecou.NSTRSS
       MFRl   = iecou.MFR1
       NCOMAT = NMATT
       jnplas = INPLAS
C*?      ifour7 = necou.IFOURB
 
       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,xmat(/1)
*       write(6,*) (xmat(iu),iu=1,xmat(/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 (jnplas.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 (jnplas.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 (MFRl.NE.1.AND.MFRl.NE.3.AND.MFRl.NE.5.AND.MFRl.NE.17.AND.
     &     MFRl.NE.31.AND.MFRl.NE.33) THEN
         KERRE = 99
         RETURN
       ENDIF
*
*  AM 5/5/00   MFRl = 33   : MODELES 19 A 24 pour le moment
*  AM 27/5/3               :  ET  44,45
       IF (MFRl.EQ.33) THEN
         IF ((jnplas.LT.19.OR.jnplas.GT.24).AND.
     &       (jnplas.NE.44.AND.jnplas.NE.45))THEN
           KERRE = 99
           RETURN
         ENDIF
       ENDIF
 
*       write(6,*) 'cconst MFR',MFRl
       IF (MFRl.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 ((jnplas.NE.29).AND.(jnplas.NE.142)) THEN
         CALL CALSIG(DEPST,DDAUX,NSTRS1,CMATE,VALMAT,VALCAR,N2EL,
     1               N2PTEL,MFRl,IFOUR7,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 (MFRl.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 ((jnplas.NE.29).AND.(jnplas.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*0.5D0
           ENDIF
10       CONTINUE
         IF (IFOUR7.EQ.-2) THEN
           SIG0(2)=0.D0
           SIG0(4)=0.D0
           IF ((jnplas.NE.29).AND.(jnplas.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 (MFRl.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 (jnplas.EQ.63.AND.MFRl.NE.1.AND.MFRl.NE.31) THEN
         KERRE=99
         RETURN
       ENDIF
       IF (DT.EQ.0.0) DT = 1.e-20
       MOTERR(1:8) = ' CONST  '
       IF (jnplas.EQ.17) THEN
         IF ((NVARI.NE.(6+4*NSTRS1)).AND.(MFRl.NE.5)) KERRE = 146
         IF ((NVARI.NE.(10+4*NSTRS1)).AND.(MFRl.EQ.5)) KERRE = 146
         IF (IFOUR7.NE.-2.AND.NCOMAT.LT.24) KERRE = 146
         IF (IFOUR7.EQ.-2.AND.NCOMAT.LT.25) KERRE = 146
         XMAX=XMAT(8)
         GOTO 30
       ENDIF
 
*       write(6,*) 'cconst inplas ifourb ncomat',jnplas,ifour7,ncomat
       IF (MFRl.NE.33) THEN
         IF (jnplas.EQ.19.AND.IFOUR7.NE.-2.AND.NCOMAT.LT. 8)KERRE = 146
         IF (jnplas.EQ.19.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT. 9)KERRE = 146
         IF (jnplas.EQ.20.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.18)KERRE = 146
         IF (jnplas.EQ.20.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.19)KERRE = 146
         IF (jnplas.EQ.21.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.12)KERRE = 146
         IF (jnplas.EQ.21.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.13)KERRE = 146
         IF (jnplas.EQ.22.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.12)KERRE = 146
         IF (jnplas.EQ.22.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.14)KERRE = 146
         IF (jnplas.EQ.23.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.16)KERRE = 146
         IF (jnplas.EQ.23.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.17)KERRE = 146
         IF (jnplas.EQ.24.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.11)KERRE = 146
         IF (jnplas.EQ.24.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.12)KERRE = 146
         IF (jnplas.EQ.25.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.26)KERRE = 146
         IF (jnplas.EQ.25.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.27)KERRE = 146
C         IF (jnplas.EQ.29.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.13)KERRE = 146
C         IF (jnplas.EQ.29.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.14)KERRE = 146
         IF (jnplas.EQ.44.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.20)KERRE = 146
         IF (jnplas.EQ.44.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.21)KERRE = 146
         IF (jnplas.EQ.45.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.27)KERRE = 146
         IF (jnplas.EQ.45.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.28)KERRE = 146
         IF (jnplas.EQ.53.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.28)KERRE = 146
         IF (jnplas.EQ.53.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.29)KERRE = 146
         IF (jnplas.EQ.61.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.18)KERRE = 146
         IF (jnplas.EQ.61.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.19)KERRE = 146
         IF (jnplas.EQ.63.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.32)KERRE = 146
         IF (jnplas.EQ.63.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.33)KERRE = 146
         IF (jnplas.EQ.70.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.14)KERRE = 146
         IF (jnplas.EQ.70.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.15)KERRE = 146
         IF (jnplas.EQ.76.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.26)KERRE = 146
         IF (jnplas.EQ.76.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.27)KERRE = 146
         IF (jnplas.EQ.77.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.18)KERRE = 146
         IF (jnplas.EQ.77.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.19)KERRE = 146
         IF (jnplas.EQ.83.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.15)KERRE = 146
         IF (jnplas.EQ.83.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.17)KERRE = 146
         IF (jnplas.EQ.84.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.13)KERRE = 146
         IF (jnplas.EQ.84.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.14)KERRE = 146
         IF (jnplas.EQ.85.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.19)KERRE = 146
         IF (jnplas.EQ.85.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.20)KERRE = 146
         IF (jnplas.EQ.86.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.17)KERRE = 146
         IF (jnplas.EQ.86.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.18)KERRE = 146
         IF (jnplas.EQ.102.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.25)KERRE = 146
         IF (jnplas.EQ.102.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.26)KERRE = 146
         IF (jnplas.EQ.130.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.10)KERRE = 146
         IF (jnplas.EQ.130.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.11)KERRE = 146
         IF (jnplas.EQ.136.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.10)KERRE = 146
         IF (jnplas.EQ.136.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.11)KERRE = 146
         IF (jnplas.EQ.137.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.10)KERRE = 146
         IF (jnplas.EQ.137.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.11)KERRE = 146
         IF (jnplas.EQ.138.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.10)KERRE = 146
         IF (jnplas.EQ.138.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.11)KERRE = 146
         IF (jnplas.EQ.139.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.10)KERRE = 146
         IF (jnplas.EQ.139.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.11)KERRE = 146
C
       ELSE
C
C        cas MFRl=33
C
         IF (jnplas.EQ.19.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.16)KERRE = 146
         IF (jnplas.EQ.19.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.17)KERRE = 146
         IF (jnplas.EQ.20.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.26)KERRE = 146
         IF (jnplas.EQ.20.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.27)KERRE = 146
         IF (jnplas.EQ.21.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.20)KERRE = 146
         IF (jnplas.EQ.21.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.21)KERRE = 146
         IF (jnplas.EQ.22.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.20)KERRE = 146
         IF (jnplas.EQ.22.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.22)KERRE = 146
         IF (jnplas.EQ.23.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.24)KERRE = 146
         IF (jnplas.EQ.23.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.25)KERRE = 146
         IF (jnplas.EQ.24.AND.IFOUR7.NE.-2.AND.NCOMAT.LT.19)KERRE = 146
         IF (jnplas.EQ.24.AND.IFOUR7.EQ.-2.AND.NCOMAT.LT.20)KERRE = 146
C
       ENDIF
C
       IF (IFOUR7.EQ.1) THEN
         KERRE = 194
         MOTERR(1:8) = 'FLUAGE'
       ENDIF
       XMAX=XMAT(5)
       IF ((jnplas.EQ.25).OR.(jnplas.EQ.53)) XMAX=XMAT(7)
       IF ((jnplas.EQ.76).OR.(jnplas.EQ.77)) XMAX=XMAT(7)
       IF ((jnplas.EQ.70).OR.(jnplas.EQ.107)) XMAX=XMAT(1)*1.D-3
       IF (jnplas.EQ.29) THEN
         CALL DERTRA(NYSMAX,YSMAX,TF,XMAX,XMAXV,TO,TO)
       ENDIF
       IF (jnplas.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 (jnplas.EQ.29) THEN
           CALL DERTRA(NYOG,YOG,TF,XMAX,XMAXV,TO,TO)
           XMAX=XMAX*1.D-3
         ELSEIF (jnplas.EQ.142) THEN
           CALL DERTRA(NYOG1,YOG1,TF,XMAX,XMAXV,TO,TO)
           XMAX=XMAX*1.D-3
         ELSEIF ((jnplas.NE.29).AND.(jnplas.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  A PARTIR DE MAINTENANT, LES DEFORMATIONS DE CISAILLEMENT 
C  NE SONT PLUS 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 (jnplas=29 ET jnplas=142).
C
C  SEULES LES FORMULATIONS SUIVANTES SONT ACCEPTEES PAR CONSTI:
C      MFRl=1  (MASSIF)
C      MFRl=5  (COQUES EPAISSES)
C      MFRl=3  (COQUES MINCES)
C      MFRl=17 (TUYAUX FISSURES)
C      MFRl=31 (BBAR)
C      MFRl=33 (POREUX)
C
       IF ((jnplas.NE.29).AND.(jnplas.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 (MFRl.EQ.1.OR.MFRl.EQ.31.OR.MFRl.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 (MFRl.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 (MFRl.EQ.3) THEN
           IF ((IFOUR7.EQ.1).OR.(IFOUR7.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 (MFRl=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 ((jnplas.NE.29).AND.(jnplas.NE.142)) THEN
           DSPT(I) = DSIGT(I)/DT
         ENDIF
40     CONTINUE
 
       IF (jnplas.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,IFOUR7,VAR0,NVARI,XMAT,NCOMAT,YUNGV0,ENUV0,
     &               XCARB,ICARA,MFRl,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,IFOUR7,VAR0,NVARI,XMAT,NCOMAT,YUNGV1,ENUV1,
     &               XCARB,ICARA,MFRl,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-TREFA)-ALFAF*(TF-TREFA)
C
         CALL ZDANUL(DSPT,NSTRS1)
         AA=CTEPS
         DO 45 I=1,NSTRS1
           IF (I.GT.3) AA=0.D0
           r_z = DEPST(I)
           DO 46 J=1,NSTRS1
             r_z = r_z + (DDINV(I,J)-DD(I,J))*SIG0(J)
46         CONTINUE
           DSPT(I)=(r_z-AA)/DT
45       CONTINUE
C
       ELSEIF (jnplas.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,IFOUR7,VAR0,NVARI,XMAT,NCOMAT,YUNGV0,ENUV0,
     &               XCARB,ICARA,MFRl,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,IFOUR7,VAR0,NVARI,XMAT,NCOMAT,YUNGV1,ENUV1,
     &               XCARB,ICARA,MFRl,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-TREFA)-ALFAF*(TF-TREFA)
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 ---------------------------------------------------------------------
       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 (MFRl.EQ.17.AND.jnplas.NE.19) GOTO 999
       IF (MFRl.NE.3.AND.MFRl.NE.17.AND.jnplas.EQ.17) THEN
         CALL XXPT1(SIG,EPSV,VAR,EVP1,VARP1,XPM1,XMAT,NSTRS1,NVARI,
     &              NCOMAT,MFRl)
       ELSE IF (MFRl.NE.3.AND.MFRl.NE.17.AND.((jnplas.GE.19.AND.
     &          jnplas.LE.24).OR.jnplas.EQ.61.OR.jnplas.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,
     &               jnplas,NCOMAT,MFRl,FII0,FII1,TI1)
       ELSE IF (MFRl.NE.3.AND.MFRl.NE.17.AND.jnplas.EQ.25) THEN
         CALL INCRE2(SIG,VAR,EVP1,VARP1,XMAT,XMAT0,NSTRS1,MFRl,NVARI,
     &               NCOMAT,VALMAT,VALMA0,DT0)
       ELSE IF ((MFRl.EQ.1.OR.MFRl.EQ.31).AND.jnplas.EQ.130) THEN
         CALL INCREP(SIG,VAR,EVP1,VARP1,XMAT,NSTRS1,MFRl,
     &               NVARI,NCOMAT)
       ELSE IF ((MFRl.EQ.1.OR.MFRl.EQ.31).AND.
     &(jnplas.EQ.136.OR.jnplas.EQ.137.OR.
     &jnplas.EQ.138.OR.jnplas.EQ.139)) THEN
         CALL INCVAR(SIG,VAR,EVP1,VARP1,XMAT,NSTRS1,MFRl,
     &               NVARI,NCOMAT,jnplas)
       ELSE IF (MFRl.NE.3.AND.MFRl.NE.17.AND.jnplas.EQ.76) THEN
         CALL INCRA2(SIG,VAR,EVP1,VARP1,XMAT,NSTRS1,MFRl,NVARI,NCOMAT)
       ELSE IF (MFRl.NE.3.AND.MFRl.NE.17.AND.jnplas.EQ.77) THEN
         CALL INCRB2(SIG,VAR,EVP1,VARP1,XMAT,NSTRS1,MFRl,NVARI,NCOMAT)
       ELSE IF (MFRl.NE.3.AND.MFRl.NE.17.AND.jnplas.EQ.53) THEN
         CALL INCRE5(SIG,VAR,EVP1,VARP1,XMAT,NSTRS1,MFRl,NVARI,NCOMAT)
       ELSE IF (MFRl.NE.3.AND.jnplas.EQ.63) THEN
         CALL INCRE7(SIG,VAR,DSPT,EVP1,VARP1,XMAT,NSTRS1,MFRl,NVARI,
     &               NCOMAT,IFOUR7)
       ELSE IF (MFRl.EQ.3.AND.((jnplas.GE.19.AND.jnplas.LE.24)
     &   .OR.jnplas.EQ.61)) THEN
         CALL INCRE3(TAU,SIG,EPSV,VAR,XMAT,EVP1,VARP1,ALFA,NSTRS1,
     &               NVARI,jnplas,NCOMAT)
       ELSE IF (MFRl.EQ.3.AND.jnplas.EQ.25) THEN
         CALL INCRE4(SIG,VAR,EVP1,VARP1,XMAT,ALFA,NSTRS1,NVARI,NCOMAT)
       ELSE IF (MFRl.EQ.3.AND.jnplas.EQ.76) THEN
         CALL INCRA4(SIG,VAR,EVP1,VARP1,XMAT,ALFA,NSTRS1,NVARI,NCOMAT)
       ELSE IF (MFRl.EQ.3.AND.jnplas.EQ.77) THEN
         CALL INCRB4(SIG,VAR,EVP1,VARP1,XMAT,ALFA,NSTRS1,NVARI,NCOMAT)
       ELSE IF (MFRl.EQ.3.AND.jnplas.EQ.53) THEN
         CALL INCRE6(SIG,VAR,EVP1,VARP1,XMAT,ALFA,NSTRS1,NVARI,NCOMAT)
       ELSE IF (MFRl.NE.3.AND.(jnplas.EQ.85.OR.jnplas.EQ.86.
     &          OR.jnplas.EQ.84.OR.jnplas.EQ.102)) THEN
         CALL DEVFLO(jnplas,SIG,EPSV,VAR,XMAT,NCOMAT,NSTRS1,NVARI,EVP1,
     &               VARP1,TAU)
       ELSE IF (MFRl.NE.3.AND.jnplas.EQ.70) THEN
         CALL INCRE8(SIG,VAR,T0,TF,EVP1,VARP1,XMAT,NSTRS1,MFRl,NVARI,
     &               NCOMAT)
C
       ELSE IF (MFRl.NE.3.AND.MFRl.NE.17.AND.jnplas.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,NCOUR7)
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,
     &                  NCOUR7)
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,
     &                  NCOUR7)
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-TREFA))+(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,MFRl,
     &                 XCARB,ICARA,IFOUR7,2,TI0,TPOINT,TINF,TSUP,ITEST,
     &                 ITHHER,TRUC,NCOUR7)
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,NCOUR7)
         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 (MFRl.NE.3.AND.MFRl.NE.17.AND.jnplas.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-TREFA))+(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,MFRl,XCARB,
     &                 ICARA,IFOUR7,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 (MFRl.EQ.17.AND.jnplas.EQ.19) THEN
         CALL TUFINC(TAU,SIG,EPSV,VAR,XMAT,XCARB,EVP1,VARP1,NSTRS1,
     &               NVARI,jnplas,NCOMAT,KERREU1)
       ELSE IF ((MFRl.EQ.1.OR.MFRl.EQ.31.OR.MFRl.EQ.33)
     &.AND.jnplas.EQ.44) THEN
         CALL POUDRA(SIG,EPSV,VAR,EVP1,VARP1,XMAT,NSTRS1,NVARI,NCOMAT,
     &               KERRE)
       ELSE IF ((MFRl.EQ.1.OR.MFRl.EQ.31.OR.MFRl.EQ.33)
     &.AND.jnplas.EQ.45) THEN
         CALL POUDRB(SIG,EPSV,VAR,EVP1,VARP1,XMAT,NSTRS1,NVARI,NCOMAT,
     &               KERRE)
 
C------ CHABOCHE
       ELSE IF ((MFRl.EQ.1.OR.MFRl.EQ.31.OR.MFRl.EQ.33)
     &          .AND.jnplas.EQ.165) THEN
         DO I=1,NCOMAT
            ZMAT(I) = XMAT(I)
         ENDDO
         CALL INCRE9(SIG,VAR,ZMAT,NSTRS1,EVP1,VARP1)
 
C------ MERIC_CAILLETAUD
       ELSE IF ((MFRl.EQ.1.OR.MFRl.EQ.31.OR.MFRl.EQ.33)
     &          .AND.jnplas.EQ.192) THEN
 
C        write(6,*) 'cconst: ncomat',ncomat
         DO I=1,NCOMAT
            ZMAT(I) = XMAT(I)
C        write(6,*) 'cconst: xmat(i)',i,xmat(i)
         ENDDO
 
* Normes vecteurs orientation
         XNN1=0.D0
         XNB1=0.D0
         IN1 = 11
         IB1 = IN1+IDIM
         DO IM=1,IDIM
           XNN1=XNN1+ZMAT(IN1+IM)**2
           XNB1=XNB1+ZMAT(IB1+IM)**2
         ENDDO
         XNN1=SQRT(XNN1)
         XNB1=SQRT(XNB1)
         IF (XNN1.LT.XZPREC.OR.XNB1.LT.XZPREC) THEN
           CALL ERREUR(239)
           RETURN
         ENDIF
 
* Definition du tenseur d'orientation cristalline
         DO K=1,6
           XM(K)=0.D0
         ENDDO
         IJ=0
         DO IM=1,IDIM
           IJ=IJ+1
           XM(IJ)=(ZMAT(IN1+IM)*ZMAT(IB1+IM))/(XNN1*XNB1)
         ENDDO
         IF (IDIM.GT.1) THEN
           DO IM=1,IDIM-1
             IJ=IJ+1
             XM(IJ)=ZMAT(IN1+IM)*ZMAT(IB1+IM+1)
             XM(IJ)=XM(IJ)+ZMAT(IN1+IM+1)*ZMAT(IB1+IM)
             XM(IJ)=0.5D0*XM(IJ)
           ENDDO
         ENDIF
C        write(6,*)'XM=',(XM(kk),kk=1,NSTRS1)
 
* Calcul increment :
         CALL INCR10(SIG,VAR,ZMAT,NSTRS1,EVP1,VARP1,XM)
 
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 (MFRl.EQ.3) GOTO 150
         IF (MFRl.EQ.17) GOTO 210
         IF (jnplas.EQ.17) THEN
           tau2=tau*0.5d0
           CALL AVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,EVP1,VARP1,
     &                XPM1,DSPT,XMAT,NSTRS1,NVARI,NCOMAT,IFOUR7,IVTEST,
     &                MFRl)
           CALL XXPT1(SIG1,EPSV1,VAR1,EVP2,VARP2,XPM2,XMAT,NSTRS1,
     &                NVARI,NCOMAT,MFRl)
           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,IFOUR7,IVTEST,
     &                MFRl)
           CALL XXPT1(SIG12,EPSV12,VAR12,EVP3,VARP3,XPM3,XMAT,NSTRS1,
     &                NVARI,NCOMAT,MFRl)
           CALL AVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,EVP3,
     &                VARP3,XPM3,DSPT,XMAT,NSTRS1,NVARI,NCOMAT,IFOUR7,
     &                IVTEST,MFRl)
           CALL XXPT1(SIG13,EPSV13,VAR13,EVP4,VARP4,XPM4,XMAT,NSTRS1,
     &                NVARI,NCOMAT,MFRl)
           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,IFOUR7,
     &                IVTEST,MFRl)
           CALL XXPT1(SIGf,EPinf,VARf,EVP4,VARP4,XPM4,XMAT,NSTRS1,
     &                NVARI,NCOMAT,MFRl)
           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,IFOUR7,IVTEST,MFRl)
C---------
         ELSE IF (jnplas.EQ.44) THEN
           CALL AVANP(TAU,SIG,EPSV,VAR,SIGf,EPinf,VARf,DSPT,EVP1,
     &                VARP1,XMAT,NSTRS1,NVARI,IFOUR7,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,IFOUR7,NCOMAT)
*           tau2=tau*0.5d0
*           CALL AVANP(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
*     &                VARP1,XMAT,NSTRS1,NVARI,IFOUR7,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,IFOUR7,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,IFOUR7,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,IFOUR7,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,IFOUR7,NCOMAT)
*C---------
         ELSE IF (jnplas.EQ.45) THEN
           CALL AVANP(TAU,SIG,EPSV,VAR,SIGf,EPinf,VARf,DSPT,EVP1,
     &                VARP1,XMAT,NSTRS1,NVARI,IFOUR7,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,IFOUR7,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,IFOUR7,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,IFOUR7,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,IFOUR7,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,IFOUR7,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,IFOUR7,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 ((jnplas.GE.19.AND.jnplas.LE.24).OR.jnplas.EQ.61.OR.
     &            jnplas.EQ.107) THEN
           tau2=tau*0.5d0
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
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 (jnplas.EQ.107) VAR1(3)=VAR(3)
           CALL INCRE1(TAU,SIG1,EPSV1,VAR1,EVP2,VARP2,XMAT,NSTRS1,
     &                 NVARI,jnplas,NCOMAT,MFRl,FII0,FII1,TI1)
           DO I=1,NSTRS1
             EVP2(I) = 0.5D0*(EVP1(I)+EVP2(I))
           ENDDO
           IF (jnplas.EQ.24.OR.jnplas.EQ.107) THEN
             IF (jnplas.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,IFOUR7,jnplas,MFRl)
           IF (jnplas.EQ.107) VAR12(3)=VAR(3)
           CALL INCRE1(TAU2,SIG12,EPSV12,VAR12,EVP3,VARP3,XMAT,NSTRS1,
     &                 NVARI,jnplas,NCOMAT,MFRl,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,IFOUR7,jnplas,MFRl)
           CALL INCRE1(TAU2,SIG13,EPSV13,VAR13,EVP4,VARP4,XMAT,NSTRS1,
     &                 NVARI,jnplas,NCOMAT,MFRl,FII12,FII1,TI1)
           DO I=1,NSTRS1
             EVP4(I) = 0.5D0*(EVP3(I)+EVP4(I))
           enddo
           IF (jnplas.EQ.24.OR.jnplas.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,IFOUR7,jnplas,MFRl)
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 (jnplas.EQ.107) VARf(3)=VAR(3)
           CALL INCRE1(TAU2,SIGf,EPinf,VARf,EVP4,VARP4,XMAT,NSTRS1,
     &                 NVARI,jnplas,NCOMAT,MFRl,FII0,FII12,TI12)
           DO I=1,NSTRS1
             EVP2(I) = (EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           enddo
           IF (jnplas.EQ.24.OR.jnplas.EQ.107) THEN
             DO I=1,NVARI
               VARP2(I)=(VARP1(I)+VARP4(I))/6.d0+VARP3(I)*2.d0/3.d0
             enddo
             IF (jnplas.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,IFOUR7,jnplas,MFRl)
           IF(jnplas.EQ.107) VARf(3)=VAR1(3)
C---------
         ELSE IF (jnplas.EQ.25) THEN
           tau2=tau*0.5d0
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           CALL INCRE2(SIG1,VAR1,EVP2,VARP2,XMAT,XMAT0,NSTRS1,MFRl,
     &                 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,IFOUR7,jnplas,MFRl)
           CALL INCRE2(SIG12,VAR12,EVP3,VARP3,XMAT,XMAT0,NSTRS1,MFRl,
     &                 NVARI,NCOMAT,VALMAT,VALMA0,DT0)
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           CALL INCRE2(SIG13,VAR13,EVP4,VARP4,XMAT,XMAT0,NSTRS1,MFRl,
     &                 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,IFOUR7,jnplas,MFRl)
           CALL INCRE2(SIGf,VARf,EVP4,VARP4,XMAT,XMAT0,NSTRS1,MFRl,
     &                 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,IFOUR7,jnplas,MFRl)
C---------
         ELSE IF (jnplas.EQ.130) THEN
           tau2=tau*0.5d0
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           CALL INCREP(SIG1,VAR1,EVP2,VARP2,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
           CALL INCREP(SIG12,VAR12,EVP3,VARP3,XMAT,NSTRS1,MFRl,NVARI,
     &                 NCOMAT)
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           CALL INCREP(SIG13,VAR13,EVP4,VARP4,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
           CALL INCREP(SIGf,VARf,EVP4,VARP4,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
C---------
C--------- LOI CHAB_SINH_R, CHAB_SINH_X (Runge Kuta 2-3)
         ELSE IF (jnplas.EQ.136.OR.jnplas.EQ.137.OR.
     a            jnplas.EQ.138.OR.jnplas.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,IFOUR7,jnplas,MFRl)
C        write(6,*)'sortie de advac'
C        write(6,*)'seconde entrée dans increp'
           CALL INCVAR(SIG1,VAR1,EVP2,VARP2,XMAT,NSTRS1,MFRl,NVARI,
     &                 NCOMAT,jnplas)
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,IFOUR7,jnplas,MFRl)
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,MFRl,NVARI,
     &                 NCOMAT,jnplas)
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,IFOUR7,jnplas,MFRl)
C        write(6,*)'sortie de advac'
C        write(6,*)'quatrième entrée dans increp'
           CALL INCVAR(SIG13,VAR13,EVP4,VARP4,XMAT,NSTRS1,MFRl,NVARI,
     &                 NCOMAT,jnplas)
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,IFOUR7,jnplas,MFRl)
C        write(6,*)'sortie de advac'
C        write(6,*)'cinquième entrée dans increp'
           CALL INCVAR(SIGf,VARf,EVP4,VARP4,XMAT,NSTRS1,MFRl,NVARI,
     &                 NCOMAT,jnplas)
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,IFOUR7,jnplas,MFRl)
C        write(6,*)'sortie de advac'
C---------
C---------
         ELSE IF (jnplas.EQ.76) THEN
           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIGf,EPinf,VARf,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           CALL INCRA2(SIGf,VARf,EVP2,VARP2,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
*           tau2=tau*0.5d0
*           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
*     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
*           CALL INCRA2(SIG1,VAR1,EVP2,VARP2,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
*           CALL INCRA2(SIG12,VAR12,EVP3,VARP3,XMAT,NSTRS1,MFRl,NVARI,
*     &                 NCOMAT)
*           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
*     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
*           CALL INCRA2(SIG13,VAR13,EVP4,VARP4,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
*           CALL INCRA2(SIGf,VARf,EVP4,VARP4,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
C---------
         ELSE IF (jnplas.EQ.77) THEN
           CALL ADVAN1(TAU,SIG,EPSV,VAR,SIGf,EPinf,VARf,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           CALL INCRB2(SIGf,VARf,EVP2,VARP2,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
 
*           tau2=tau*0.5d0
*           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
*     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
*           CALL INCRB2(SIG1,VAR1,EVP2,VARP2,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
*           CALL INCRB2(SIG12,VAR12,EVP3,VARP3,XMAT,NSTRS1,MFRl,NVARI,
*     &                 NCOMAT)
*           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
*     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
*           CALL INCRB2(SIG13,VAR13,EVP4,VARP4,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
*           CALL INCRB2(SIGf,VARf,EVP4,VARP4,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
C---------
         ELSE IF (jnplas.EQ.53) THEN
           tau2=tau*0.5d0
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           CALL INCRE5(SIG1,VAR1,EVP2,VARP2,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
           CALL INCRE5(SIG12,VAR12,EVP3,VARP3,XMAT,NSTRS1,MFRl,NVARI,
     &                 NCOMAT)
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           CALL INCRE5(SIG13,VAR13,EVP4,VARP4,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
           CALL INCRE5(SIGf,VARf,EVP4,VARP4,XMAT,NSTRS1,MFRl,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,IFOUR7,jnplas,MFRl)
C---------
         ELSE IF (jnplas.EQ.85.OR.jnplas.EQ.86.OR.jnplas.EQ.84.OR.
     &       jnplas.EQ.102) THEN
           tau2=tau*0.5d0
           CALL ADVFLO(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,NCOMAT,
     &                 MFRl)
           CALL DEVFLO(jnplas,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,IFOUR7,jnplas,NCOMAT,
     &                 MFRl)
           CALL DEVFLO(jnplas,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,IFOUR7,jnplas,
     &                 NCOMAT,MFRl)
           CALL DEVFLO(jnplas,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,IFOUR7,jnplas,
     &                 NCOMAT,MFRl)
           CALL DEVFLO(jnplas,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,IFOUR7,jnplas,NCOMAT,
     &                 MFRl)
C----------
         ELSE IF (jnplas.EQ.63) THEN
           tau2=tau*0.5d0
           CALL ADVDDI(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,NCOMAT,
     &                 MFRl)
           CALL INCRE7(SIG1,VAR1,DSPT,EVP2,VARP2,XMAT,NSTRS1,MFRl,
     &                 NVARI,NCOMAT,IFOUR7)
           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,IFOUR7,jnplas,NCOMAT,
     &                 MFRl)
           CALL INCRE7(SIG12,VAR12,DSPT,EVP3,VARP3,XMAT,NSTRS1,MFRl,
     &                 NVARI,NCOMAT,IFOUR7)
           CALL ADVDDI(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,
     &                 NCOMAT,MFRl)
           CALL INCRE7(SIG13,VAR13,DSPT,EVP4,VARP4,XMAT,NSTRS1,MFRl,
     &                 NVARI,NCOMAT,IFOUR7)
           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,IFOUR7,jnplas,
     &                 NCOMAT,MFRl)
           CALL INCRE7(SIGf,VARf,DSPT,EVP4,VARP4,XMAT,NSTRS1,MFRl,
     &                 NVARI,NCOMAT,IFOUR7)
           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,IFOUR7,jnplas,NCOMAT,
     &                 MFRl)
C--------
         ELSE IF (jnplas.EQ.70) THEN
           tau2=tau*0.5d0
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           t12=(t0+tf)*0.5d0
           CALL INCRE8(SIG1,VAR1,T0,T12,EVP2,VARP2,XMAT,NSTRS1,MFRl,
     &                 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,IFOUR7,jnplas,MFRl)
           CALL INCRE8(SIG12,VAR12,T12,TF,EVP3,VARP3,XMAT,NSTRS1,MFRl,
     &                 NVARI,NCOMAT)
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           CALL INCRE8(SIG13,VAR13,T12,TF,EVP4,VARP4,XMAT,NSTRS1,MFRl,
     &                 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,IFOUR7,jnplas,MFRl)
           CALL INCRE8(SIGf,VARf,T12,TF,EVP4,VARP4,XMAT,NSTRS1,MFRl,
     &                 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,IFOUR7,jnplas,MFRl)
C---------
         ELSE IF (jnplas.EQ.29) THEN
C
 143       TAU2=0.5D0*TAU
*           write(6,*) 'tau=',tau
*           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)
*         write(6,*) 'varp1(3)=',varp1(3)
*       Do i=1,6
*         write(6,*) 'sig1(',i,')=',sig1(i)
*       enddo
*       Do i=1,6
*         write(6,*) 'epsv1(',i,')=',epsv1(i)
*       enddo
*       Do i=1,7
*         write(6,*) 'var1(',i,')=',var1(i)
*       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,NCOUR7)
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,
     &                  NCOUR7)
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,
     &                  NCOUR7)
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-TREFA))+(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,
     &                 MFRl,XCARB,ICARA,IFOUR7,2,TI12,TPOINT,TINF,TSUP,
     &                 ITEST,ITHHER,TRUC,NCOUR7)
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,
     &                 MFRl,XCARB,ICARA,IFOUR7,2,TI12,TPOINT,TINF,TSUP,
     &                 ITEST,ITHHER,TRUC,NCOUR7)
           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,NCOUR7)
           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,
     &                  NCOUR7)
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,
     &                  NCOUR7)
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-TREFA))+(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,
     &                 MFRl,XCARB,ICARA,IFOUR7,2,TI1,TPOINT,TINF,TSUP,
     &                 ITEST,ITHHER,TRUC,NCOUR7)
           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,
     &                 MFRl,XCARB,ICARA,IFOUR7,2,TI1,TPOINT,TINF,TSUP,
     &                 ITEST,ITHHER,TRUC,NCOUR7)
           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)
*       Do i=1,6
*         write(6,*) 'sig1(',i,')=',sig1(i)
*       enddo
*       Do i=1,6
*         write(6,*) 'epsv1(',i,')=',epsv1(i)
*       enddo
*       Do i=1,3
*         write(6,*) 'var1(',i,')=',var1(i)
*       enddo
c           IF (TAU.NE.T) GOTO 143
C---------
         ELSE IF (jnplas.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-TREFA))+(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,MFRl,
     &                 XCARB,ICARA,IFOUR7,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,MFRl,
     &                 XCARB,ICARA,IFOUR7,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-TREFA))+(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,MFRl,
     &                 XCARB,ICARA,IFOUR7,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,MFRl,
     &                 XCARB,ICARA,IFOUR7,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 (jnplas.EQ.165) THEN
           tau2=tau*0.5d0
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,ZMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           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,IFOUR7,jnplas,MFRl)
           CALL INCRE9(SIG12,VAR12,ZMAT,NSTRS1,EVP3,VARP3)
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,ZMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           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,IFOUR7,jnplas,MFRl)
           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,IFOUR7,jnplas,MFRl)
 
C--------- MERIC_CAILLETAUD
         ELSE IF (jnplas.EQ.192) THEN
           tau2=tau*0.5d0
           CALL ADVAN1(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,ZMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           DO I=1,NCOMAT
             ZMAT(I)  = 0.5D0*(XMAT(I)+XMAT0(I))
           ENDDO
           CALL INCR10(SIG1,VAR1,ZMAT,NSTRS1,EVP2,VARP2,XM)
           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,IFOUR7,jnplas,MFRl)
           CALL INCR10(SIG12,VAR12,ZMAT,NSTRS1,EVP3,VARP3,XM)
           CALL ADVAN1(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,ZMAT,NSTRS1,NVARI,IFOUR7,jnplas,MFRl)
           CALL INCR10(SIG13,VAR13,ZMAT,NSTRS1,EVP4,VARP4,XM)
           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,IFOUR7,jnplas,MFRl)
           DO I=1,NCOMAT
             ZMAT(I)  = XMAT(I)
           ENDDO
           CALL INCR10(SIGF,VARF,ZMAT,NSTRS1,EVP4,VARP4,XM)
           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,IFOUR7,jnplas,MFRl)
C---------
         ENDIF
         GOTO 250
C _____________________________________________________________________
150    CONTINUE
C----------------------------------------------------------------------
C  CALCULATIONS FOR GENERALISED STRESS/STRAIN FORMULATIONS
C----------------------------------------------------------------------
         IF ((jnplas.GE.19.AND.jnplas.LE.24).OR.jnplas.EQ.61) THEN
           tau2=tau*0.5d0
           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,NCOMAT)
           CALL INCRE3(TAU2,SIG1,EPSV1,VAR1,XMAT,EVP2,VARP2,ALFA,
     &                 NSTRS1,NVARI,jnplas,NCOMAT)
           DO I=1,NSTRS1,1
             EVP2(I) = 0.5d0*(EVP1(I)+EVP2(I))
           ENDDO
           IF (jnplas.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,IFOUR7,jnplas,NCOMAT)
           CALL INCRE3(TAU2,SIG12,EPSV12,VAR12,XMAT,EVP3,VARP3,ALFA,
     &                 NSTRS1,NVARI,jnplas,NCOMAT)
           CALL ADVAN2(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &                 EVP3,VARP3,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,
     &                 NCOMAT)
           CALL INCRE3(TAU2,SIG13,EPSV13,VAR13,XMAT,EVP4,VARP4,ALFA,
     &                 NSTRS1,NVARI,jnplas,NCOMAT)
           DO I=1,NSTRS1
              EVP4(I) = 0.5d0*(EVP3(I)+EVP4(I))
           enddo
           IF (jnplas.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,IFOUR7,jnplas,
     &                 NCOMAT)
           CALL INCRE3(TAU2,SIGf,EPinf,VARf,XMAT,EVP4,VARP4,ALFA,
     &                 NSTRS1,NVARI,jnplas,NCOMAT)
           DO I=1,NSTRS1,1
             EVP2(I) =(EVP1(I)+EVP4(I))/6.d0+EVP3(I)*2.d0/3.d0
           enddo
           IF (jnplas.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,IFOUR7,jnplas,NCOMAT)
C---------
         ELSE IF (jnplas.EQ.25) THEN
           tau2=tau*0.5d0
           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,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,IFOUR7,jnplas,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,IFOUR7,jnplas,
     &                 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,IFOUR7,jnplas,
     &                 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,IFOUR7,jnplas,NCOMAT)
C---------
         ELSE IF (jnplas.EQ.76) THEN
           CALL ADVAN2(TAU,SIG,EPSV,VAR,SIGf,EPinf,VARf,DSPT,
     &                 EVP1,VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,
     &                 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,IFOUR7,jnplas,NCOMAT)
*           tau2=tau*0.5d0
*           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,
*     &                 EVP1,VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,
*     &                 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,IFOUR7,jnplas,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,IFOUR7,jnplas,
*     &                 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,IFOUR7,jnplas,
*     &                 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,IFOUR7,jnplas,NCOMAT)
C---------
         ELSE IF (jnplas.EQ.77) THEN
           CALL ADVAN2(TAU,SIG,EPSV,VAR,SIGf,EPinf,VARf,DSPT,EVP1,
     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,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,IFOUR7,jnplas,NCOMAT)
*           tau2=tau*0.5d0
*           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,
*     &                 VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,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,IFOUR7,jnplas,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,IFOUR7,jnplas,
*     &                 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,IFOUR7,jnplas,
*     &                 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,IFOUR7,jnplas,NCOMAT)
C---------
         ELSE IF (jnplas.EQ.53) THEN
           tau2=tau*0.5d0
           CALL ADVAN2(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,
     &                 EVP1,VARP1,XMAT,NSTRS1,NVARI,IFOUR7,jnplas,
     &                 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,IFOUR7,jnplas,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,IFOUR7,jnplas,
     &                 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,IFOUR7,jnplas,
     &                 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,IFOUR7,jnplas,NCOMAT)
C---------
         ENDIF
         GOTO 250
C
C CAS D'UN TUYAU FISSURE EN FLUAGE SUIVANT LA LOI "NORTON"
C
 210     CONTINUE
         IF (jnplas.NE.19) GOTO 999
         tau2=tau*0.5d0
         CALL TUFADV(TAU2,SIG,EPSV,VAR,SIG1,EPSV1,VAR1,DSPT,EVP1,VARP1,
     &               XMAT,XCARB,NSTRS1,NVARI,IFOUR7,jnplas,NCOMAT)
         CALL TUFINC(TAU2,SIG1,EPSV1,VAR1,XMAT,XCARB,EVP2,VARP2,NSTRS1,
     &               NVARI,jnplas,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,IFOUR7,jnplas,
     &               NCOMAT)
         CALL TUFINC(TAU2,SIG12,EPSV12,VAR12,XMAT,XCARB,EVP3,VARP3,
     &               NSTRS1,NVARI,jnplas,NCOMAT,KERREU1)
         CALL TUFADV(TAU2,SIG12,EPSV12,VAR12,SIG13,EPsv13,VAR13,DSPT,
     &               EVP3,VARP3,XMAT,XCARB,NSTRS1,NVARI,IFOUR7,jnplas,
     &               NCOMAT)
         CALL TUFINC(TAU2,SIG13,EPSV13,VAR13,XMAT,XCARB,EVP4,VARP4,
     &               NSTRS1,NVARI,jnplas,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,IFOUR7,jnplas,
     &               NCOMAT)
         CALL TUFINC(TAU2,SIGf,EPinf,VARf,XMAT,XCARB,EVP4,VARP4,NSTRS1,
     &               NVARI,jnplas,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,IFOUR7,jnplas,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 (MFRl.EQ.17 .AND. jnplas.EQ.19) THEN
         RA=SQRT(XX(1)**2 + XX(6)**2)/ERRABS
       ENDIF
       SQRA = SQRT(RA)
C
       IF (jnplas.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 (jnplas.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 (jnplas.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 (jnplas.EQ.165.OR.jnplas.EQ.192) THEN
C        Convergence variables internes Chaboche ou Meric_Cailletaud
         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 ((jnplas.EQ.29).OR.(jnplas.EQ.142)) TAU =MIN(TAU,TAUX)
           DELTAT=TPOINT*TAU
           TI1=TI0+DELTAT
           GOTO 80
         ELSEIF ( RA.GT.(BORNE)) THEN
           TAU = TAU/SQRA
           IF ((jnplas.EQ.29).OR.(jnplas.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 (jnplas.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,NCOUR7)
 
         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 (jnplas.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 ((jnplas.EQ.29).OR.(jnplas.EQ.107).OR.(jnplas.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 (jnplas.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 (MFRl.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 (jnplas=29).
C
C  SEULES LES FORMULATIONS SUIVANTES SONT ACCEPTEES PAR CONSTI:
C      MFRl=1  (MASSIF)
C      MFRl=5  (COQUES EPAISSES)
C      MFRl=3  (COQUES MINCES)
C      MFRl=17 (TUYAUX FISSURES)
C      MFRl=31 (BBAR)
C      MFRl=33 (POREUX)
C
       IF ((jnplas.NE.29).AND.(jnplas.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 (MFRl.EQ.1.OR.MFRl.EQ.31.OR.MFRl.EQ.33) THEN
           DO 14 I=4,NSTRS1
             EPIN0(I)=EPIN0(I)*2.D0
             EPINF(I)=EPINF(I)*2.D0
 14        CONTINUE
C
C  Cas des coques épaisses
C  Les termes de cisaillement apparaissent
C  au delà de la deuxieme composante
C
         ELSE IF (MFRl.EQ.5) THEN
           DO 15 I=3,NSTRS1
             EPIN0(I)=EPIN0(I)*2.D0
             EPINF(I)=EPINF(I)*2.D0
 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 (MFRl.EQ.3) THEN
           IF ((IFOUR7.EQ.1).OR.(IFOUR7.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 (MFRl=17)
C
         ENDIF
       ENDIF
C
C===========================================================
 998  CONTINUE
      RETURN
C
C     On detecte des NaN ou des Inf dans les contraintes
 259  CONTINUE
      CALL ERREUR(460)
 
      RETURN
      END