Sources Esope | Cast3M

ccoinc
C CCOINC    SOURCE    OF166741  25/11/04    21:15:09     12349          
      SUBROUTINE CCOINC(wrk52,wrk53,wrk54,wrk2,wrk3,
     &  IB,IGAU,NBPGAU,ecou,necou,iecou)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
C---------------------------------------------------------------------
C               ECOULEMENT PLASTIQUE POUR UN POINT
C                       D 'APRES INCA
C---------------------------------------------------------------------
C
C      EN ENTREE :
C
C                  SIG0       CONTRAINTES AU DEBUT DU PAS
C                  DEPST      INCREMENT DE DEFORMATIONS TOTALES
C                               ( THERMIQUE ENLEVEE )
C                  var0     VARIABLES INTERNES DEDUT DU PAS
C                  bid     VARIABLES EXTERNES DEBUT DU PAS
C                  VAREXF     VARIABLES EXTERNES FIN DU PAS
C                  XMAT       COEFFICIENTS DU MATERIAU
C                  PRECIS     PRECISION POUR ITERATIONS INTERNES
C                  WORK      DES CARACTERISTIQUES
C                  TRAC       COURBE DE TRACTION
C                  MFR1        INDICE DE FORMULATION
C                  NSTRS1      NOMBRE DE CONTRAINTES CA2000
C                  INPLAS     NUMERO DU MODELE DE PLASTICITE
*                  DDAUX = MATRICE DE HOOKE ELASTIQUE
*                  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
*                  IFOU  = OPTION DE CALCUL
*                  IB    = NUMERO DE L ELEMENT COURANT
*                  IGAU  = NUMERO DU POINT COURANT
*                  EPAIST= EPAISSEUR
*                  NBPGAU= NBRE DE POINTS DE GAUSS
*                  MELE  = NUMERO DE L ELEMENT FINI
*                  NPINT = NBRE DE POINTS D INTEGRATION
*                  NBGMAT= NBRE DE POINTS DANS SEGMENT DE CARACTERISTIQUES
*                  NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES
*                  SECT  = SECTION
*                  LHOOK = TAILLE DE LA MATRICE DE HOOKE
*                  TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI = TABLEAUX UTILISES
*                  UTILISES POUR LE CALCUL DE LA MATRICE DE HOOKE
C
C     EN SORTIE :
C
C                  SIGF        CONTRAINTES A LA FIN DU PAS
C                  varf      VARIABLES INTERNES FIN DU PAS
C                  DEFP        INCR. DE DEFORMATIONS PLASTIQUES
C                  KERRE       CODE D'ERREUR
C                   = 0  SI TOUT OK
C                   = 99 SI FORMULATION NON DISPONIBLE
C    EN PLASTICITE
C                   = 1  SI DEPS EST NEGATIF
C                   = 2  SI NOMBRE MAX D'ITERATIONS INTERNES DEPASSE
C    EN FLUAGE
C                   = 3  DIMINUER LE PAS ( FLUAGE RUNGE )
C                   = 4  DEPS EST NEGATIF
C                   = 5  PAS DE CONVERGENCE INTERNE ( ANCIEN NSTOP=0 )
C                   = 6  MATRICE SINGULIERE DANS CHAFLU
C    EN PLUS
C                   = 7  PROBLEME DANS LES CARACT DU TUYAU
C                   = 49 PB DANS LES ITER. INTERNES DRUCKER-PRAGER
C                   = 51 PB DANS LES ITER. INTERNES DRUCKER-PRAGER
C                   = 52 PB DANS LES ITER. INTERNES DRUCKER-PRAGER
C                   = 53 PB DANS LES ITER. INTERNES DRUCKER-PRAGER
C                   = 75 PB AVEC LA COURBE DE TRACTION
C
C-----------------------------------------------------------------------
C      VARIABLES PASSEES PAR LES COMMON ECOU ET NECOU
C
C  IFOUR     INDICE DU TYPE DE PROBLEME
C            -2  CONTRAINTES PLANES
C            -1  DEFORMATIONS PLANES
C             0
C             N  SERIE DE FOURIER
C  IMAPLA    INDICE DE MATERIAU PLASTIQUE
C            0  MATERIAU ELASTIQUE
C            4  MODELE DE CHABOCHE
C            5  MODELE DE DRUCKER-PRAGER
C            7  MODELE DE ROUSSELIER
C  ITYP      TYPE DE FORMULATION MECANIQUE
C       ITYP=1   CAS DES ELEMENTS MASSIFS
C       ITYP=2   CAS DES COQUES
C       ITYP=3   CAS DES MEMBRANES
C       ITYP=4   CAS DES CABLES ET DES BARRES
C       ITYP=5   CAS QUELCONQUE
C       ITYP=6   CAS DES CONTRAINTES PLANES
C       ITYP=7   CAS DES COQUES A NU=0. OU CONTRAINTES PLANES
C       ITYP=8   CAS DES MEMBRANES A NU=0. OU CONTRAINTES PLANES                 C       ITYP=9   CAS DES COQUES EPAISSES
C       ITYP=10  CAS DES JOINTS
C       ITYP=11  CAS DES POUTRES
C       ITYP=12  CAS DES TUYAUX
C       ITYP=13  CAS DES COQUES AVEC CISAILLEMENT TRANSVERSE
C  IBI       1 SI ERREUR DANS TRACTI , 0 SINON
C  IPLAST    1 SI ON A PLASTIFIE , 0 SINON
C  ICINE  = 1 ECROUISSAGE CINEMATIQUE  ( CHABOCHE COMPRIS )
C  JFLUAG = 1 ON FLUE AVEC SIGMA
C  JFLUAG = 2 ON FLUE AVEC (SIG-X)
C  LFLUAG = 0 ON ECROUIT EN CAS DE FLUAGE
C  LFLUAG = 1 ON N'ECROUIT PAS EN CAS DE FLUAGE
C  JNTRIN = 0 LE TEMPS INTRINSEQUE EST LE TEMPS REEL
C  JNTRIN = 1 ON EVOLUE LE TEMPS INTRINSEQUE EN PLASTICITE ET FLUAGE
C  JNTRIN = 2 ON N'EVOLUE LE TEMPS INTRINSEQUE QU'EN FLUAGE
C  JNTRIN = 3 ON EST EN PLASTICITE ISOCHRONE
C
C  TETI = TEMPERATURE AU DEBUT DU PAS
C  TET  = TEMPERATURE A LA FIN DU PAS
C-----------------------------------------------------------------------
 
-INC PPARAM
-INC CCOPTIO
-INC DECHE
 
-INC TECOU
 
      SEGMENT WRK2
        REAL*8 TRAC(LTRAC)
      ENDSEGMENT
 
      SEGMENT WRK3
        REAL*8 WORK(LW),WORK2(LW2)
      ENDSEGMENT
 
      REAL*8   DDOT
      EXTERNAL DDOT
 
      DIMENSION SIGMA(6),DSIGMA(6),SPHER(6),AUXIL(6),COEF(40)
      DIMENSION DIV(7),SIGFIN(6),DEFPLA(6)
      DIMENSION TABID(6),TABID2(6)
      REAL*8 CRIGI(12)
      REAL*8 EI(1)
C ZZZZZZZZZZZZZZZZZZZZZZZ
C    DIMENSIONS A REVOIR
C ZZZZZZZZZZZZZZZZZZZZZ
      DIMENSION DHOOK(1),SIG(130),EPS(130)
      DIMENSION ORMAT(1),ANORM(1)
      DIMENSION DSIGZE(6),DETFNO(1),NNNAA(13)
      DATA ITMAX/15/
      DATA PI4,R33,R22/0.785398164D0,1.732050808D0,1.414213562D0/
      DATA A,B,C,D/.577350269D0,.7071067814D0,.4082482904D0,
     . -0.8164965808D0/
      DATA A1/1.D0/
      DATA A2/.5D0/
      DATA A3/3.D0/
      DATA COEF/0.D0,1.5D0,1.5D0,3.D0,3.D0,3.D0,0.5D0,1.5D0,3.D0,0.5D0,
     . 1.5D0,3.D0,0.D0,0.D0,1.D0,0.D0,0.D0,0.D0,1.D0,
     . 1.D0,1.D0,2.D0,2.D0,2.D0,0.5D0,1.5D0,0.D0,3.D0,
     . 1.D0,0.D0,0.D0,1.D0,1.D0,1.D0,0.5D0,1.5D0,0.D0,3.D0,3.D0,3.D0/
      DATA XRF/0.0001D0/
      DATA NNNAA/6,6,3,3,6,4,6,1,6,3,6,6,6/
C-----------------------------------------------------------------------
C         CONVENTION DE REMPLISSAGE DES MEMOIRES
C-----------------------------------------------------------------------
C
C     XMAT(1)    MODULE D'YOUNG
C     XMAT(2)    COEFFICIENT DE POISSON
C
C     TRAC(1 A 2*NCOURB) COURBE DE TRACTION
C     WORK(  " +1) ALFAH POUR COQUES PLASTICITE GLOBALE
C     WORK(  " +..) DONNEES POUR CRITERE POUTRES GLOBALES
C
C    MODELE ISOTROPE
C            var0(1)  EPS*
C
C    MODELE CINEMATIQUE LINEAIRE
C            var0(1)  EPS*
C            var0(2) A var0(1+IBOU)  DEPLACEMENT DE LA SPHERE
C
C    MODELE CHABOCHE
C            XMAT(5) ....      COEFFICIENTS A,C,...
C            var0(1)  EPS*
C            var0(2) A var0(1+IBOU) DEPLACEMENT DE LA SPHERE 1
C            var0(2+IBOU) A var0(1+2*IBOU)  "    "  "   "    2
C            ICENT2 = 1 SI ON A 2 CENTRES
C
C    MODELE DRUCKER-PRAGER
C            XMAT(5) .... XMAT(13) = LES CONSTANTES DU MODELE
C            var0(1)  EPS*
C
C    MODELE ROUSSELIER
C
C    MODELE DE CYCLAGE
C
C    MODELE A TEMPS INTRINSEQUE ( FLUAGE ET PLASTICITE ISOCHRONE )
C            var0(.)  TIMEXI   TEMPS INTRINSEQUE DE DEBUT DE PAS
C
C--------------------------------------------------------------------
C                     REMPLISSAGE
C-----------------------------------------------------------------------
      YUNG=XMAT(1)
      XNU =XMAT(2)
C
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
C    REMPLISSAGE A ACTUALISER SELON LES MODELES
CZZZZZZZZZZZZZZZZZZZZZZZZZZZ
CZZZZZZZZZZZZZZZZZZZZZZ
C      VARIABLES A RECUPERER
CZZZZZZZZZZZZZZZZZZZZZZ
C  ORMAT  MATERIAU ORTHOTROPE      DANS XMAT
C  ANORM    A VOIR         DANS XMAT OU DANS VARIN
C  EPENT    DANS VARIN
C  VECPRO,VALPRO  DANS XMAT
C  COVNMS,CVNMSD  DANS XMAT ( SELON LES CAS )
C
C  IT    DANS VAREXT
C  EPSFLU DANS VARINT
C  SIG,EPS DANS XMAT
C   DETFNO  DANS VARIN ????
C  DSIGZE DANS VARIN  ( POUR FLUAGE RELAX )
C  SIPLAD DANS VARIN  ( POUR FLUAGE RELAX  )  ...
C  ... ( CE SONT LES DEFORMATIONS OU CONTRAINTES (A VOIR) PLASTIQUES
C                            AU DEBUT DU PAS )
C  TEMPS , HPAS DANS VAREXT
C  TET,TETI DANS VAREXT
C
      EPSM1=var0(1)
      TIMEXI=var0(1)
      DPSM1=var0(1)
      DPSM2=var0(1)
      EPENT=var0(1)
      EPSFLU=var0(1)
      DETFNO(1)=var0(1)
      DSIGZE(1)=var0(1)
      TEMPS0=bid(1)
      IT=nint(bid(3))
      TEMPS=bid(1)
      TETI=bid(2)
      TET =bid(2)
      DO I=1,NCOURB
        SIG(I)=TRAC(2*I-1)
        EPS(I)=TRAC(2*I)
      ENDDO
      ORMAT(1)=XMAT(1)
      DHOOK(1)=XMAT(1)
C----------------------------------------------------------------------
C                       INITIALISATIONS
C----------------------------------------------------------------------
      KERRE=0
      JA=1
      JC=1
      IA=1
 
      IFOURL = necou.IFOURB
      MFRL   = iecou.MFR1
      NSTRS1 = iecou.NSTRSS
      jNPLAS = INPLAS
C
C  PETIT TEST SUR NU POUR CERTAINS CAS
C
      IF(MFRL.EQ.3.AND.IFOURL.EQ.-2.AND.XNU.NE.0.D0) THEN
        KERRE=38
        RETURN
      ENDIF
CZZZZZZZZZZZZZ
C  PROVISOIRE
CZZZZZZZZZZZZ
      ANORM(1)=XMAT(5)
      ICENT2=0
      IF (jNPLAS.EQ.12.OR.jNPLAS.EQ.13) ICENT2=1
      IF (jNPLAS.EQ.7) NUMCHA=1
      IF (jNPLAS.EQ.11) NUMCHA=2
      IF (jNPLAS.EQ.12) NUMCHA=3
      IF (jNPLAS.EQ.13) NUMCHA=4
C
      ECTEST=0.5D00*PRECIS
      HPAS=TEMPS-TEMPS0
      MCOD=1
 
C cas particulier pour elastique non-lineaire equiplas
 
        if (jNPLAS.EQ.87) THEN
             EPST=0.D0
             SIG0(1)=0.D0
             DEPST(1)=DEPST(1) + VAR0(2)
             EPSM1 =0.D0
       ENDIF
*
*   CALCUL DE DSIGT
*
      CALL CALSIG(DEPST,DDAUX,NSTRS1,CMATE,VALMAT,VALCAR,
     1        N2EL,N2PTEL,MFRL,IFOURL,IB,IGAU,EPAIST,
     2      NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,
     3    XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD)
*
      IF(IRTD.NE.1) THEN
         KERRE=69
         GOTO 1000
      ENDIF
*
*    TRAITEMENT DES COQ2 ET MASSIF 2D EN CAS DE MATERIAU
*    UNIDIRECTIONNEL
*
*      LES VARIABLES INTERNES SONT DANS LE REPERE UNIDIRECTIONNEL
*      ON SE LIMITE AU CAS AXISYMETRIQUE ?
*
            IF (CMATE.EQ.'UNIDIREC'.AND.jNPLAS.NE.0) THEN
              IF(MELE.EQ.44.OR.(MFRL.EQ.1.AND.IFOURL.LE.0)) THEN
                 APHI=0.D0
*
*                   CAS MASSIF : ON AJOUTE L'ANGLE DU REPERE
*                   DE L'ELEMENT AVEC LE REPERE GLOBAL
*
                 IF(MFRL.EQ.1) APHI=ATAN2(TXR(2,1),TXR(1,1))
                 DO 1996 I=1,NSTRS1
                    TABID(I)=SIG0(I)
                    TABID2(I)=DSIGT(I)
 1996            CONTINUE
                 ANG = ATAN2(XMAT(4),XMAT(3))  + APHI
                 CALL CHREP2(TABID,SIG0,ANG,NSTRS1,MFRL,1,KERRE)
                 CALL CHREP2(TABID2,DSIGT,ANG,NSTRS1,MFRL,1,KERRE)
                 IF(KERRE.NE.0) RETURN
               ELSEIF(MFRL.EQ.1.AND.IFOURL.EQ.2) THEN
                 DO 996 I=1,NSTRS1
                    TABID(I)=SIG0(I)
                    TABID2(I)=DSIGT(I)
 996             CONTINUE
                 CALL CHREP3(TABID,SIG0,XMAT,TXR,NSTRS1,1,1,KERRE)
                 CALL CHREP3(TABID2,DSIGT,XMAT,TXR,NSTRS1,1,1,KERRE)
 
              ENDIF
            ENDIF
*
      CALL VISAVI(SIG0,DSIGT,var0,SIGMA,DSIGMA,SPHER,AUXIL,
     .            SIGF,DEFP,varf,SIGFIN,DEFPLA,
     .   DSIGZE,ICENT2,MCOD,IBOU,MFRL,NSTRS1,WORK,CMATE,ecou,necou)
      IF(ITYP.EQ.0) THEN
      KERRE=99
      RETURN
      ENDIF
C-----------------------------------------------------------------------
      EI(1)=YUNG
      UNSYG=1.D00/YUNG
      U5SNU=1.5D00/(1.D00+XNU)
      U5DMU=U5SNU*YUNG
      SREF=YUNG*XRF
      IPLAST=0
      LCGDF=36
      LPLUS=6
      ITER=1
      LAPOIN=0
C
C  CAS DES COQUES EN GLOBAL - ON RECUPERE LE ALFAH
C
      ALFAH=1.D0
      IF(ITYP.EQ.2) ALFAH=WORK(2)**2
      IF(ITYP.EQ.7) ALFAH=WORK(2)**2
      UNALF=0.D0
      IF(ALFAH.GE.1.D-12) UNALF=1.D0/ALFAH
C
C  CAS DES BARRES ET DES CABLES
C
      IF(ITYP.EQ.4) THEN
      DIV(1)=1.D0
      DIV(2)=1.D0
      DIV(3)=1.D0/WORK(1)
      ENDIF
 
C
C  CAS DES POUTRES
C
      IF(ITYP.NE.11) GO TO 841
      DIV(1)=1.D0/WORK(4)
      DIV(2)=1.D0
      DIV(3)=1.D0
      DIV(4)=WORK(10)/WORK(1)
      DIV(5)=WORK(11)/WORK(2)
      DIV(6)=WORK(12)/WORK(3)
      IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(WORK(1)*WORK(4))
      IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(WORK(2)*WORK(4))
      IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(WORK(3)*WORK(4))
      GO TO 761
 841  CONTINUE
C
C  CAS DES TUYAUX
C
      IF(ITYP.NE.12) GO TO 842
      EPAIS=WORK(1)
      REXT=WORK(2)
      RMOY=REXT-EPAIS*0.5D0
      RACO=WORK(3)
      GAM=1.D0
      XK=1.D0
      IF(RACO.EQ.0.D0) GO TO 765
      XLAM=RMOY*RMOY/EPAIS/RACO
      XK=1.65D0*XLAM
      IF(XK.LT.1.D0) XK=1.D0
      GAM=0.8888888888888889D0*(XLAM)**0.6666666666666667D0
      IF(GAM.LT.1.D0) GAM=1.D0
 765  CONTINUE
C
C     NB 23/09/98
C     VALEURS PAR DEFAUT POUR LES CFFX CFMX CFMY
C     CFMZ CFPR ( COEFFICIENTS POUR CALCULER LES
C     CONTRAINTES DE MEMBRANE, TORSION, FLEXIONS
C     DANS LE PLAN, HORS PLAN ET CIRCONFERENTIELLE
C     DUE A LA PRESSION )
C     POUR L'INSTANT PAS DE CONTRAINTE CIRCONFERENTIELLE
C     DUE A LA PRESSION ON N'UTILISE DONC PAS DIV(7)
C
C
      DIV(1)=1.D0
      DIV(2)=1.D0
      DIV(3)=1.D0
      DIV(4)=R33
      DIV(5)=PI4*GAM
      DIV(6)=DIV(5)
      DIV(7)=R22*0.5D0
      IF(IDIM.EQ.2) THEN
        PRES1=WORK(6)
        CISA1=WORK(7)
        IXCAR1=12
        IDEB1=8
      ELSE IF(IDIM.EQ.3) THEN
        PRES1=WORK(7)
        CISA1=WORK(8)
        IXCAR1=13
        IDEB1=9
      ENDIF
C
      JDIV1=2
      DO IBA=IDEB1,IXCAR1
        JDIV1=JDIV1+1
        VCAR1=WORK(IBA)
        IF (VCAR1.NE.-1.D0) DIV(JDIV1)=WORK(IBA)
      ENDDO
C
C     NB 23/09/98
C     TRANSFERT DE CFFX DANS DIV(1) ET REMISE A
C     1.D0 DE DIV(3)
C
      DIV(1) = DIV(3)
      DIV(3)=1.D0
C
      VX=WORK(4)
      VY=WORK(5)
      VZ=WORK(6)
      CALL TUYCAR(WORK,CISA1,VX,VY,VZ,KERRE,2)
      IF (KERRE.NE.0) THEN
        KERRE=7
        RETURN
      ENDIF
C
C     NB 23/09/98
C     DIVISION DE DIV(5) ET DIV(6) PAR XK
C     CAR TUYCAR SORT LES INERTIES MODIFIEES
C     POUR LE CALCUL DES DEFORMATIONS PAR EPSIG
C
      DIV(5)=DIV(5)/XK
      DIV(6)=DIV(6)/XK
C
      DIV(1)=DIV(1)/WORK(4)
      DIV(4)=DIV(4)*RMOY/WORK(1)
      DIV(5)=DIV(5)*RMOY/WORK(2)
      DIV(6)=DIV(6)*RMOY/WORK(3)
761   CONTINUE
      E(1)=YUNG
      E(2)=0.D0
      E(3)=0.D0
      E(4)=YUNG*1.5D0/(1.D0+XNU)
      E(5)=YUNG
      E(6)=YUNG
  842 CONTINUE
C
      IF(ICINE.EQ.0.OR.JFLUAG.EQ.1) GO TO 204
C
C       ON EST EN CINEMATIQUE ( PLASTIQUE OU FLUAGE )
C       DANS CE CAS ON FLUE AVEC SIGMA-ALPHA QUE L'ON CALCULE
C
      DO IBA=1,IBOU
        SIGMA(IBA)=SIGMA(IBA)-SPHER(IBA)
      ENDDO
 204  CONTINUE
C
      IF(JNTRIN.EQ.0) GO TO 830
      TIMEXF=TIMEXI
      GO TO 831
830   TIMEXI=TEMPS-HPAS
      TIMEXF=TEMPS
831   CONTINUE
      IF(ITYP.NE.11.AND.ITYP.NE.12.AND.ITYP.NE.4) GO TO 844
      DO IBA=1,IBOU
      SIGMA(IBA)=SIGMA(IBA)*DIV(IBA)
      IF(ICINE.EQ.0) GO TO 845
      IF(JFLUAG.EQ.1.AND.LFLUAG.EQ.1) GO TO 845
      SPHER(IBA)=SPHER(IBA)*DIV(IBA)
845   CONTINUE
      DSIGMA(IBA)=DSIGMA(IBA)*DIV(IBA)
      ENDDO
844   CONTINUE
      DO IBA=1,IBOU
        STOT(IBA)=SIGMA(IBA)+DSIGMA(IBA)
      ENDDO
C----------------------------------------------------------------------
C               CALCUL DE LA LIMITE ELASTIQUE
C----------------------------------------------------------------------
      IF(LFLUAG.EQ.1) GO TO 261
      IF(IMAPLA.NE.4) GO TO 262
      BPSTAR=EPSM1
      ICOD=1
      CALL CHALIM(BPSTAR,SELAS,XMAT,TET,ICOD,
     . BID,BID,BID,BID,BID1,BID2,BI3,BI4,BI5,BI6,IBID,IBID,NUMCHA)
      GO TO 261
  262 IF(IMAPLA.NE.5) GO TO 263
*                     REMPLACEMENT  XMAT(7) PAR XMAT(12) AM  8/3/90
*                     TEST SELON EPSM1                   AM 28/1/91
      IF(EPSM1.EQ.0.D0) THEN
      SELAS = XMAT(12)
      ELSE
        SELAS = XMAT(7) + EPSM1 * XMAT(13)
*
*                     AM  24/5/93  TEST SUR SELAS
*
        IF( SELAS.LT.0.D0) THEN
          SELAS = 0.D0
        ENDIF
*
      ENDIF
      GO TO 261
  263 CONTINUE
      EPSTAR=EPSM1
      IF(ICINE.EQ.1) EPSTAR=0.D0
      CALL TRACTI(SELAS,EPSTAR,SIG,EPS,NCOURB,2,IBI)
      IF(IBI.EQ.1) THEN
        KERRE=75
        GOTO 1000
      ENDIF
*
261   CONTINUE
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
C   A VOIR    CE QU'IL Y A DANS ANORM
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
      IF(IMAPLA.EQ.7) COVNMS(1)=ANORM(2*JC-1)
      IF(IMAPLA.EQ.7.AND.IT.EQ.1) ANORM(2*JC)=1.D-20
*
*      SI MATERIAU DRUCKER PRAGER ON CHERCHE LE CRITERE
*               AVEC LEQUEL ON DOIT FAIRE LA PROJECTION ET LE CRITERE
      IXMAT=5
      IF(IMAPLA.EQ.5.AND.EPSM1.EQ.0.D0) IXMAT=10
      CALL CKRIT(IMAPLA,STOT,ITYP,XMAT(IXMAT),ALFAH,COVNMS,XINV,SSTAR)
      SJ1=XINV(1)
      SJ2=XINV(2)
      IF(IMAPLA.EQ.0) GO TO 3
      IF(JFLUAG.NE.0) GO TO 4
C---------------------------------------------------------------------
C          A T'ON DEPASSE LA LIMITE ELASTIQUE
C---------------------------------------------------------------------
      PETI=1.1D00*ECTEST*SSTAR
      CALL EPSPRE(SSTAR,SELAS,PETI,ITRY)
      IF(ITRY.EQ.0.AND.SSTAR.GT.SELAS) GO TO 4
C----------------------------------------------------------------------
C         ON N'A PAS PLASTIFIE
C----------------------------------------------------------------------
    3 CONTINUE
      IF(ITYP.EQ.11.OR.ITYP.EQ.12.OR.ITYP.EQ.4) THEN
      DO IBA=1,IBOU
        IF(ICINE.NE.0)  SPHER(IBA)=SPHER(IBA)/DIV(IBA)
        STOT(IBA)=STOT(IBA)/DIV(IBA)
      ENDDO
      ENDIF
      DO IBA=1,IBOU
        SIGFIN(IBA)=STOT(IBA)
        DEFPLA(IBA)=0.D0
        DSIGMA(IBA)=0.D0
      ENDDO
      IF(jNPLAS.EQ.87) THEN
            VARF(2)=DEPST(1)
      ENDIF
      STARF=SSTAR
      EPST=EPSM1
      GO TO 31
C----------------------------------------------------------------------
C                ON A PLASTIFIE
C----------------------------------------------------------------------
 4    CONTINUE
      IPLAST=1
      PENTE=0.D0
      SLOPE=0.D0
      DO IBA=1,IBOU
        W1(IBA)=SIGMA(IBA)
      ENDDO
      CALL CKRIT(IMAPLA,W1,ITYP,XMAT(5),ALFAH,COVNMS,XINV,S0)
      EPSTAR=EPSM1
C
C    EN FLUAGE SI ON LA VITESSE DE FLUAGE NE DEPEND QUE
C    DU EPSILON EQUIVALENT DE FLUAGE CUMULE ON VA LE LIRE DANS EPSFLU
C
      IF(JFLUAG.NE.0.AND.IFLUPL.EQ.1) EPSTAR=EPSFLU
      EPSIN=EPSTAR
      SEL=0.D0
      SPLA=1.D0
      IF(JFLUAG.EQ.0) THEN
        CALL SIGELP(SIGMA,DSIGMA,SIGEL,DSIGP,STOT,PETI,
     .       ITYP,SEL,SPLA,IBOU,S0,SELAS,XMAT(5),COVNMS,ALFAH,
     .       IMAPLA,SSTAR)
        CALL SHIFTD(SIGEL,DALPHA,IBOU)
      ELSE
        DO IBA=1,IBOU
          SIGEL(IBA)=SIGMA(IBA)+SEL*DSIGMA(IBA)
          DALPHA(IBA)=SIGEL(IBA)
          DSIGP(IBA)=SPLA*DSIGMA(IBA)
        ENDDO
      ENDIF
      DTIMEX=SPLA*HPAS
      IF(JNTRIN.EQ.1.OR.JNTRIN.EQ.3) TIMEXF=TIMEXI+DTIMEX
C
C           CALCUL DES PENTES
C
      IF(ABS(SPLA).LT.0.95D00) IELPLA=0
      IF(ICYCL.NE.1) GO TO 71
C
C            CAS DES MODELES DE CYCLAGE
C
      CALL NORMAL(A1,A2,A3,SIGEL,W1,ALFAH,SELAS,IBOU,ITYP)
CZZZZZZZZZZZZZZZZZZZZZZ
C    ON RECUPERE LA NORMALE DE LA DERNIERE SORTIE
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
      DO IBA=1,IBOU
        W2(IBA)=ANORM(IBA)
      ENDDO
C
C   CALCUL DU PRODUIT SCALAIRE DE LA NORMALE AU POINT OU ON SORT
C   AVEC LA DERNIERE NORMALE OU ON EST SORTI
C
      X=DDOT(IBOU,W1,1,W2,1)
C                          ON STOCKE LA NORMALE
      DO IBA=1,IBOU
        ANORM(IBA)=W1(IBA)
      ENDDO
CZZZZZZZZZZZZZZZZZZZZZZZ
C   IL FAUDRA LA STOCKER DANS varf AU LIEU DE ANORM
CZZZZZZZZZZZZZZZZZZZZZZZZZ
C
C           CALCUL DU EPSI AVEC QUOI CALCULER LA PENTE
C
      EPSEN=0.5D00*(X+1.D00)*EPENT
      EPSIN=EPSEN
   71 CONTINUE
      IF(JFLUAG.EQ.0) GO TO 40
      IF(IMAPLA.EQ.4) GO TO 5
C
C         INITIALISATIONS POUR LE FLUAGE
C         QUEL EST LE EPSILON CUMULE DE FLUAGE
C
      SLOPE=0.D0
      PENTE=0.D0
      IF(ICINE.EQ.0) GO TO 5
      IF(JFLUAG.EQ.2) GO TO 40
C
C      CAS CINE + FLUAGE : CALCUL DE LA PENTE A LA COURBE DE TRACTION
C
      IF(ICINE.EQ.2) EPSIN=EPS(NCOURB)
      CALL TRACTI(PENTE,EPSIN,SIG,EPS,NCOURB,1,IBI)
      IF(IBI.EQ.1) THEN
        KERRE=75
        GOTO 1000
      ENDIF
*
      GO TO 5
   40 CONTINUE
      IF(IMAPLA.EQ.5.OR.IMAPLA.EQ.4) GO TO 5
C
C CALCUL DE LA PENTE DE LA COURBE DE TRACTION AU POINT X-EPSTAR
C
      SLOPE=0.D0
      CALL TRACTI(PENTE,EPSIN,SIG,EPS,NCOURB,1,IBI)
      IF(IBI.EQ.1) THEN
        KERRE=75
        GOTO 1000
      ENDIF
*
      IF(ICINE.EQ.0) GO TO 5
 6    CONTINUE
      SLOPE=PENTE
      PENTE=0.
      PANTIN=SLOPE
      IF(ICINE.EQ.1) GO TO 5
      EPSIN=EPS(NCOURB)
      CALL TRACTI(PANTIN,EPSIN,SIG,EPS,NCOURB,1,IBI)
      IF(IBI.EQ.1) THEN
        KERRE=75
        GOTO 1000
      ENDIF
*
      PENTE=SLOPE-PANTIN
      SLOPE=PANTIN
 5    CONTINUE
      IF(IMAPLA.NE.4) GO TO 5621
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C                       MODELE CHABOCHE
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      IF(JFLUAG.EQ.0)
     .  CALL CHABOK(YUNG,XNU,IA,EI,
     .  XMAT,SPHER,JA,IBOU,SI,DEPS,EPST,EPSTAR,WORK2,KERRE,SN,
     .  AUXIL,NUMCHA,ecou,necou)
      IF(JFLUAG.NE.0)
     .  CALL CHAFLU(YUNG,XNU,IA,EI,SSTAR,
     1  XMAT,SPHER,IBOU,SI,DEPS,EPST,EPSTAR,WORK2,
     2  AUXIL,DPSM1,DPSM2,KERRE,NUMCHA,ecou,necou)
      IF(KERRE.NE.0) RETURN
      IF(JFLUAG.NE.0.AND.IMAPLA.EQ.4) GO TO 689
      IF(IMAPLA.EQ.4) GO TO 570
 5621 CONTINUE
      IF(IMAPLA.NE.5) GO TO 9997
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C                MODELE DRUCKER-PRAGER  -  PRELIMINAIRES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     ON EST DANS LE CAS DU MATERIAU DRUCKER-PRAGER.ON VA FAIRE UN
C     TEST POUR SAVOIR SI L'ON DOIT PROJETER AU SOMMET DU CONE ET
C     EVENTUELLEMENT FAIRE CETTE PROJECTION.
C     DANS CETTE VERSION,LE TEST CONSISTE A SAVOIR SI L'ETAT DE
C     CONTRAINTES CALCULE ELASTIQUEMENT SE TROUVE D'UN COTE OU DE
C     L'AUTRE D'UN PLAN ORTHOGONAL A L'AXE ET PASSANT PAR LE SOMMET
C     DU CONE,PUIS A SAVOIR SI L'ON EST A L'INTERIEUR OU A
C     L'EXTERIEUR DU CONE DES NORMALES AU CRITERE DE DRUCKER.
C     ON CALCULE LE I1 QUI DELIMITE LES DEUX REGIONS ET ON LE
C     COMPARE AVEC LE PREMIER INVARIANT DU TENSEUR CALCULE ELASTIQUEMENT
C
C
C     SNN=XMAT(7)+XMAT(13)*EPSTAR
*
*                     AM  24/5/93  TEST SUR SNN
*
*     IF( SNN*XMAT(7).LT.0.D0) THEN
*       SNN=0.D0
*       XMAT(13)=0.D0
*     ENDIF
*
C
C     MILL  17/3/93
C     SI XMAT(5)=0. PAS DE TEST PAR RAPPORT AU SOMMET
C
C     IF(XMAT(5).EQ.0.D0) GO TO 99970
C
C     XI1LIM=SNN/XMAT(5)
C     IF(SJ1.LE.XI1LIM) GO TO 99970
C
C     ON TESTE ENSUITE PAR RAPPORT AU CONE DES NORMALES
C
C     TESTDR=-XMAT(9)*SJ1/(2.D0*XMAT(8))+SJ2+XI1LIM*XMAT(9)/
C    .                    (2.D0*XMAT(8)*XMAT(5))
C     IF(TESTDR.GT.0.D00)GO TO 99970
C
C     ON EST MAINTENANT DANS LE CAS OU L'ON DOIT PROJETER AU SOMMET
C
C     CALL SOMDRU(IBOU,SI,DEPS,EPST,EPSTAR,SNN,XMAT(5),YUNG,XNU,KERRE)
C     IF(KERRE.NE.0) RETURN
C     LAPOIN=1
C     GO TO 570
99970 CONTINUE
      CALL SHIFTD(DHOOK,WORK2(109),36)
      CALL PRJDRU(SSTAR,SELAS,PENTE,IBOU,SI,DEPS,EPST,EPSTAR,ITER,
     .   SN,WORK2(109),WORK2,WORK2(37),WORK2(73),XMAT(5),YUNG,
     .   XNU,LAPOIN,KERRE,ecou,necou)
      IF(KERRE.NE.0) RETURN
      GO TO 570
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C                        MODELE ROUSSELIER
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9997  IF(IMAPLA.EQ.7) THEN
      CALL FISSUR
C     CALL FISSUR(YUNG,XNU,U5DMU,XMAT(5),A,B,C,ANORM,
C    1 JC,SIG,EPS,SELAS,SSTAR,IBOU,JKI,PENTE,EPSTAR,EPST,ITMAX,COEF,
C    1 IA,S0,SI,DEPS,SN,IMPRI)
      GO TO 570
      ENDIF
C
      IF (IRELAX.NE.0) THEN
C
C  MODIF DU 18/3/86  ON MET SIPLAD A 0. ( AVANT ON AVAIT DEFP )
C
        DO IBA=1,IBOU
          SIPLAD(IBA)=0.D00
          DSIGP0(IBA)=DSIGZE(IBA)
        ENDDO
      ENDIF
      IF(JGRDEF.LE.1) GO TO 420
C
C       INTEGRATION EN GRANDES DEFORMATIONS D N'EST PLUS ISOTROPE
C ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
C            A VOIR|||||||||||||||||||||||||
      CALL SHIFTD(DHOOK,WORK2,36)
      RAI=DETFNO(1)
      CALL T6CDCT
C     CALL T6CDCT(CGDFNO((IA-1)*LCGDF+1),WORK2,4,W1,RAI)
      CALL PRJGDF
C     CALL PRJGDF(SSTAR,SELAS,PENTE,IBOU,SI,DEPS,EPST,EPSTAR,
C    1 ITER,SN,TET,WORK2,WORK2(37),WORK2(73),PANTIN,TIMEXF,SIG,EPS,IA,
C    1 WORK2(109))
      CALL SHIFTD(W1,DEFPLA,6)
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
      GO TO 570
  420 CONTINUE
C----------------------------------------------------------------------
C                 ECRITURE SOUS FORME DIAGONALE
C               DIAGONALISATION DE SIGEL ET DSIGP
C----------------------------------------------------------------------
      NK=1
      CALL XYZDIA(DSIGP,SIGEL,W1,W2,A,B,C,D,ITYP,VECPRO)
      IF(IRELAX.NE.0)
     1 CALL XYZDIA(SIPLAD,DSIGP0,W1,W2,A,B,C,D,ITYP,VECPRO)
C
      GO TO (60,61,61,63,64,65,66,66,68,69,78,79,80),ITYP
C
   80 CONTINUE
C                 CONTRAINTES PLANES AVEC TOUS LES CISAILLEMENTS
C
      RAI=YUNG/(1.D0-XNU*XNU)
      E(1)=RAI*(1.D0+XNU)*0.5D0
      E(2)=RAI*(1.D0-XNU)*1.5D0
      E(3)=0.D0
      E(4)=E(2)
      E(5)=E(2)
      E(6)=E(2)
      JKI=34
      GO TO 615
C
  79  CONTINUE
C                   TUYAUX
  78  CONTINUE
C                   POUTRES
C     LE REMPLISSAGE DE E EST FAIT EN AMONT
C
      JKI=28
      GO TO 615
   69 CONTINUE
      GO TO 640
   68 CONTINUE
      GO TO 640
   66 CONTINUE
C                   COQUE OU MEMBRANE A NU =0 OU CP
C
      RAI=YUNG
      E(1)=RAI
      IF(ICINE.NE.0) E(1)=E(1)+SLOPE
      E(2)=0.D00
      E(3)=0.D00
      E(4)=E(1)*ALFAH
      E(5)=0.D00
      E(6)=0.D00
      W1(1)=DSIGP(1)/(E(1)+PENTE)
      IF(ITYP.EQ.7) THEN
        IF(ALFAH.NE.0.D0) THEN
        W1(4)=DSIGP(4)/(E(4)+PENTE)
        ELSE
        W1(4)=0.D0
        ENDIF
      ENDIF
      JKI=18
      DEPS=VNMISD(W1,ITYP,ALFAH,COVNMS)
      GO TO 640
   65 CONTINUE
C                      CONTRAINTES  PLANES
C
      RAI=YUNG/(1.D0-XNU*XNU)
      E(1)=RAI*(1.D0+XNU)*0.5D0
      E(2)=RAI*(1.D0-XNU)*1.5D0
      E(3)=0.D0
      E(4)=E(2)
      E(5)=E(3)
      E(6)=E(3)
      JKI=24
      GO TO 615
   64 CONTINUE
C                   MATER QUELCONQUE
C
      JKI=18
      NK=1
      DO IBA=1,6
        E(IBA)=VALPRO(IBA)
      ENDDO
      IF(ICINE.NE.0) THEN
        DO IBA=1,6
          E(IBA)=E(IBA)+SLOPE
        ENDDO
      ENDIF
      DO IBA=1,6
        r_z = E(IBA)+PENTE
        IF (r_z.NE.0.D00) THEN
          W1(IBA)=DSIGP(IBA)/r_z
        ENDIF
      ENDDO
      DEPS=VNMISD(W1,ITYP,ALFAH,CVNMSD)
      GO TO 640
C
   61 CONTINUE
      RAI=YUNG/(1.D0-XNU*XNU)
      E(1)=RAI*(1.D0+XNU)*0.5D0
      E(2)=RAI*(1.D0-XNU)*1.5D0
      E(3)=E(2)
      E(4)=E(1)*ALFAH
      E(5)=E(2)*ALFAH
      E(6)=E(3)*ALFAH
      JKI=6
  615 CONTINUE
      IF(ICINE.NE.0) THEN
        DO IK=1,6
          E(IK)=E(IK)+SLOPE
        ENDDO
      ENDIF
      IBO=2
      IF(ITYP.EQ.3.OR.ITYP.EQ.6) IBO=1
      IF(ITYP.NE.11.AND.ITYP.NE.12.AND.ITYP.NE.4) W1(3)=0.D0
      DO I=1,IBOU
        IF(ITYP.EQ.6.AND.I.EQ.3) GO TO 614
        W1(I)=0.D0
        IF(PENTE+E(I).EQ.0.) GO TO 614
        W1(I)=DSIGP(I)/(E(I)+PENTE)
  614   CONTINUE
      ENDDO
      DEPS=VNMISD(W1,ITYP,ALFAH,COVNMS)
      GO TO 640
C
   63 E(3)=U5DMU/U5SNU
      IF(ICINE.NE.0) E(3)=E(3)+SLOPE
      NK=3
      JKI=12
      E(1)=0.D0
      E(2)=0.D0
      E(4)=0.D0
      E(5)=0.D0
      E(6)=0.D0
      DEPS=DSIGP(3)/(E(3)+PENTE)
      if(imapla.eq.1) deps = abs(deps)
      GO TO 640
C
   60 CONTINUE
      JKI=0
      E(1)=0.
      DO I=2,6
        E(I)=U5DMU
      ENDDO
      IF(ICINE.NE.0) THEN
        DO I=1,6
          E(I)=E(I)+SLOPE
        ENDDO
      ENDIF
      EMU=PENTE+U5DMU
      RAI=YUNG
      IF(IMAPLA.EQ.5) EMU=U5DMU+3*XMAT(5)*XMAT(5)*RAI/(1.D0-2.D0*XNU)
      DEPS=(SSTAR-SELAS)/EMU
C
  640 CONTINUE
C
                  IF(JFLUAG.NE.0) GO TO 41
C-----------------------------------------------------------------------
C       PLASTICITE  - PROJECTION SUR LA SURFACE DE CHARGE
C-----------------------------------------------------------------------
      CALL PROJVM(SIG,EPS,EPST,EPSTAR,DEPS,PENTE,SN,PANTIN,
     1 SELAS,ITER,SO,SSTAR,SI,NK,IBOU,COEF,JKI,KERRE,ecou,necou)
C
      GO TO 57
C-----------------------------------------------------------------------
C                       FLUAGE
C-----------------------------------------------------------------------
   41 CONTINUE
      CALL FLUAGE(ITER,EPSTAR,S0,EPST,SSTAR,SI,DPSM1,DPSM2,
     . IBOU,DEPS,ICONV,TIMEXI,TIMEXF,YUNG,KERRE,ecou,necou )
      IF(KERRE.EQ.0) GO TO 8693
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
C   TESTER  KERRE EN SORTIE DE FLUAGE
C  ON RECUPERE AUUSI ICONV EN SORTIE
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
8693  CONTINUE
689   CONTINUE
C-----------------------------------------------------------------------
C          ON STOKE LE EPS DE FLUAGE DANS EPSFLU DANS TOUS LES CAS
C-----------------------------------------------------------------------
      EPSFLU=EPST
C           ON CALCULE LE EPSILON INELASTIQUE CUMULE
C
      EPST=EPSM1
      IF(LFLUAG.EQ.1) GO TO 414
C             ON ECROUIT LE RAYON DE LA SURFACE DE CHARGE
C
      EPST=DEPS+EPSM1
  414 CONTINUE
      IF(IMAPLA.EQ.4) GO TO 570
      IF(ICINE.EQ.0) GO TO 413
      IF(LFLUAG.EQ.1) GO TO 57
C            CAS CINEMATIQUE
C
      SLIM=SELAS
      GO TO 57
  413 CONTINUE
C         CAS ISOTROPE
C
      IF(LFLUAG.EQ.1) GO TO 57
C            ON CHERCHE LE NOUVEAU RAYON DE LA SPHERE
C
      CALL TRACTI(SF,EPST,SIG,EPS,NCOURB,2,IBI)
      IF(IBI.EQ.1) THEN
        KERRE=75
        GOTO 1000
      ENDIF
*
      SLIM=SF
C
  57   CONTINUE
C
C----------------------------------------------------------------------
C                         RETOUR AUX AXES X Y Z
C----------------------------------------------------------------------
      CALL DIAXYZ(ITYP,SIGT,SIGEL,VECPRO,IBO,A,B,C,B1,B2)
      IF(JFLUAG.EQ.0) GO TO 10
      IF((IT.EQ.1.AND.IRELAX.EQ.0).OR.(ICONV.NE.0.AND.IRELAX.NE.0))
     . CALL DIAXYZ(ITYP,DSIGP,STOT,VECPRO,IBO,A,B,C,B1,B2)
10    CONTINUE
C-----------------------------------------------------------------------
C                     TRAITEMENT FINAL
C-----------------------------------------------------------------------
  570 CONTINUE
      EPSM1F=EPST
      IF(ICYCL.EQ.1) EPENTF=EPSEN+DEPS
      IF(IMAPLA.NE.7) GO TO 9100
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
C     VOIR SI ON A BESOIN DE SIGT
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
      DO IBA=1,IBOU
         SIGT(IBA)=SIGT(IBA)-W1(IBA)
      ENDDO
 9100 CONTINUE
      IF(ICINE.EQ.0) GO TO 913
      IF(LFLUAG.EQ.1) GO TO 913
      IF(JFLUAG.EQ.1) GO TO 94
      IF(IMAPLA.EQ.4) GO TO 90
C----------------------------------------------------------------------
C         CAS CINEMATIQUE - CALCUL DE DALPHA
C----------------------------------------------------------------------
      Z=PANTIN*DEPS/SELAS
      IF(Z.LT.0.01D00) GO TO 910
      EZZ=1.D0-EXP(-Z)
      Z=1.D0-EZZ/Z
      EZZ=EZZ/(1.D0-Z)
      Z=Z/(1.D0-Z)
      DO IBA=1,IBOU
        DALPHA(IBA)=DALPHA(IBA)*EZZ+Z*(SIGEL(IBA)-DALPHA(IBA))
      ENDDO
      GO TO 90
  910 CONTINUE
      DO IBA=1,IBOU
        DALPHA(IBA)=Z*0.5D0*(DALPHA(IBA)+SIGEL(IBA))
      ENDDO
      GO TO 90
   94 CONTINUE
      IF(IMAPLA.EQ.4) GO TO 913
      CALL MOVCEN(PENTE,DEPS,SELAS,JA,IBOU,DALPHA,SIGMA,SPHER,SIGEL)
  90  CONTINUE
 913  CONTINUE
C----------------------------------------------------------------------
C         CAS CINEMATIQUE    COORDONNEES DU NOUVEAU CENTRE DE LA SPHERE
C----------------------------------------------------------------------
      IF(ICINE.EQ.0) GO TO 29
C
C  MODELE CHABOC  MISE A JOUR DE SPHER DANS CHABOC
C  SPHER CONTIENT LA SOMME DES DEPLACEMENTS DES SPHERES
C
       IF(IMAPLA.EQ.4) GO TO 29
      IF(LFLUAG.EQ.1.AND.JFLUAG.NE.0) GO TO 29
      DO IBA=1,IBOU
        SPHER(IBA)=SPHER(IBA)+DALPHA(IBA)
      ENDDO
   29 CONTINUE
C-----------------------------------------------------------------------
C  CALCUL DU SIGMA FINAL ET DE LA NOUVELLE LIMITE ELASTIQUE
C-----------------------------------------------------------------------
      IF(ITYP.NE.11.AND.ITYP.NE.12.AND.ITYP.NE.4) GO TO 846
      DO IBA=1,IBOU
      STOT(IBA)=STOT(IBA)/DIV(IBA)
      IF(ICINE.EQ.0) GO TO 847
      IF(LFLUAG.EQ.1.AND.JFLUAG.NE.0) GO TO 847
      SPHER(IBA)=SPHER(IBA)/DIV(IBA)
      DALPHA(IBA)=DALPHA(IBA)/DIV(IBA)
847   CONTINUE
      SIGEL(IBA)=SIGEL(IBA)/DIV(IBA)
      ENDDO
846   CONTINUE
      DO IBA=1,IBOU
        SIGFIN(IBA)=SIGEL(IBA)
      ENDDO
      STARF=SI
      IF(JFLUAG.EQ.0) SLIMF=SN
      IF(JNTRIN.NE.0) EPENTF=TIMEXF
C
      IF(DEPS.LT.0.D00) KERRE=1
C          CALCUL DU NOUVEAU SIGMAP
      DO IBA=1,IBOU
        DSIGMA(IBA)=STOT(IBA)-SIGEL(IBA)
      ENDDO
      IF(ICINE.EQ.0) GO TO 240
      IF(JFLUAG.EQ.1) GO TO 240
      IF(LFLUAG.EQ.1) GO TO 240
      DO IBA=1,IBOU
        DSIGMA(IBA)=DSIGMA(IBA)-DALPHA(IBA)
      ENDDO
  240 CONTINUE
*
      IF(ITYP.EQ.6.AND.IMAPLA.EQ.5.AND.LAPOIN.EQ.0) THEN
      DO IBA=1,IBOU
        DEFPLA(IBA)=EPSPLA(IBA)
      ENDDO
      EPST=EPSM1 + VONEPS(DEFPLA,ITYP,ALFAH,COVNMS)
      GO TO 31
      ENDIF
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
C    ATTENTION LE EPSIG NE MARCHE QUE POUR CERTAINS CAS
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
      CALL EPSIG(DSIGMA,DEFPLA,IFOURL,YUNG,XNU,ITYP,ORMAT,WORK)
      IF(ITYP.EQ.6.AND.IMAPLA.NE.5) DEFPLA(3)=-DEFPLA(1)-DEFPLA(2)
  31  CONTINUE
C-----------------------------------------------------------------------
C                  REMPLISSAGE
C-----------------------------------------------------------------------
      varf(1)=EPST
      if(jNPLAS.eq.87) THEN
          DO I=1,IBOU
              VARF(1+I)=DEPST(I)
         ENDDO
      endif
CZZZZZZZZZZZ
C PROVISOIRE
CZZZZZZZZZZZ
      IF(JFLUAG.EQ.1) varf(1)=EPSFLU
      IF(ICINE.EQ.0) GO TO 280
C              ON CALCULE LA VRAIE CONTRAINTE EN CINEMATIQUE
      DO 281 IBA=1,IBOU
      IF(JFLUAG.EQ.1) GO TO 284
      SIGFIN(IBA)=SIGFIN(IBA)+SPHER(IBA)
284   IF(LFLUAG.EQ.1) GO TO 281
      IF(ICENT2.NE.0) SPHER(IBA)=SPHER(IBA)-AUXIL(IBA)
  281 CONTINUE
  280 CONTINUE
*
*  MISE A 0. DES COMPOSANTES NON CALCULEES SI ITYP=8
*
      IF(ITYP.EQ.8.AND.IBOU.LT.NSTRS1) THEN
         DO 285 IBA=2,NSTRS1
            DEFPLA(IBA)=0.D0
            SIGFIN(IBA)=0.D0
              IF(ICENT2.NE.0) SPHER(IBA)=0.D0
 285     CONTINUE
      ENDIF
*
      MCOD=2
      CALL VISAVI(SIG0,DSIGT,var0,SIGMA,DSIGMA,SPHER,AUXIL,
     .            SIGF,DEFP,varf,SIGFIN,DEFPLA,
     .   DSIGZE,ICENT2,MCOD,IBOU,MFRL,NSTRS1,WORK,CMATE,ecou,necou)
C
C      RETRAITEMENT POUR LES COQ2 ET MASSIFS  EN MATERIAU
C      UNIDIRECTIONNEL
C
 
              IF (CMATE.EQ.'UNIDIREC'.AND.jNPLAS.NE.0) THEN
                 IF(MELE.EQ.44.OR.(MFRL.EQ.1.AND.IFOURL.LE.0)) THEN
                  DO 1997 I=1,NSTRS1
                    TABID(I)=SIGF(I)
                    TABID2(I)=DEFP(I)
 1997             CONTINUE
                  CALL CHREP2(TABID,SIGF,-ANG,NSTRS1,MFRL,1,KERRE)
                  CALL CHREP2(TABID2,DEFP,-ANG,NSTRS1,MFRL,2,KERRE)
                  IF(KERRE.NE.0) RETURN
                 ELSEIF(MFRL.EQ.1.AND.IFOURL.EQ.2) THEN
                  DO 997 I=1,NSTRS1
                    TABID(I)=SIGF(I)
                    TABID2(I)=DEFP(I)
 997              CONTINUE
                  CALL CHREP3(TABID,SIGF,XMAT,TXR,NSTRS1,1,2,KERRE)
                  CALL CHREP3(TABID2,DEFP,XMAT,TXR,NSTRS1,2,2,KERRE)
                 ENDIF
              ENDIF
*
 1000 RETURN
      END