C VISCOE    SOURCE    PV        21/01/15    21:15:01     10850          
      SUBROUTINE VISCOE(CM,EPSIPP,EPSE,DTT,DCT,DT
     &                      ,DELTAT,DTTDT,DCTDT,DTDT)
C
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      DIMENSION EPSIPP(3),EPSILT(3),CM(*),DDHOOK(6,6),
     1 SIGP(3),SIGPT(3),SIGPC(3)
      PARAMETER (XUnDemi=0.5D0,XZERO=0.D0,UN=1.D0,DEUX=2.D0)
      PARAMETER (XPETIT=1.D-12)

C
C     recuperation des donnees du modele
C

        YOUN=CM(1)
        XNU=CM(2)
        XMDt=CM(14)
        XNDt=CM(15)
        At=CM(16)
        Bt=CM(17)
        XMDc=CM(18)
        XNDc=CM(19)
        Ac=CM(20)
        Bc=CM(21)
        ED0=CM(22)
C
C     Matrice de hooke (pour le 3D)
C
        AUX0=YOUN/((Un+XNU)*(Un-XNU-XNU))
        AUX=AUX0*(Un-XNU)
        AUX1=AUX0*XNU
        GEGE=XUnDemi*YOUN/(Un+XNU)
        DDHOOK(1,1)=AUX
        DDHOOK(2,1)=AUX1
        DDHOOK(3,1)=AUX1
        DDHOOK(1,2)=AUX1
        DDHOOK(2,2)=AUX
        DDHOOK(3,2)=AUX1
        DDHOOK(1,3)=AUX1
        DDHOOK(2,3)=AUX1
        DDHOOK(3,3)=AUX
        DDHOOK(4,4)=GEGE
        DDHOOK(5,5)=GEGE
        DDHOOK(6,6)=GEGE
C
C
C      on calcule l'endommagement et les contraintes
C
C
C        on calcule ALFAT ALFAC DT et DC puis D
C
C
        DO 200 ISTRS=1,3
           SIGP(ISTRS)= XZERO
           DO 210 JSTRS=1,3

             SIGP(ISTRS)=SIGP(ISTRS)+DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS)
210        CONTINUE
200     CONTINUE
C
C        on calcule le signe des contraintes elastiques
C
        DO 300 ISTRS=1,3
         IF (SIGP(ISTRS).LT. XZERO) THEN
            SIGPC(ISTRS)=SIGP(ISTRS)
            SIGPT(ISTRS)=XZERO
         ELSE
            SIGPT(ISTRS)=SIGP(ISTRS)
            SIGPC(ISTRS)=XZERO
         END IF
300     CONTINUE
        TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3)
        TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)
C
C        on calcule les deformations dues aux contraintes positives
C
        DO 400 ISTRS=1,3
          EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN
400     CONTINUE
C
C        on en deduit ALFAT et ALFAC
C
        ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) +
     1          MAX(XZERO,EPSIPP(2))*EPSILT(2) +
     2          MAX(XZERO,EPSIPP(3))*EPSILT(3)
        ALFAT = ALFAT/(EPSE*EPSE)
        ALFAC = UN - ALFAT
C
C       modification pour la bi ou tricompression
C
        IF (TRSIGC.LT. -XPETIT .AND. TRSIGT.LT.XPETIT) THEN
          GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+SIGPC(3)*SIGPC(3)
          GAMMA=-SQRT(GAMMA)/TRSIGC
          EPSE=EPSE*GAMMA
        END IF
C
C        on calcule DT et DC puis D
**      dpot=((EPSE-ED0-(UN/At)*(dtt/(UN-dtt))**(UN/Bt))/XMdt)**XNdt
        expt=(dtt/(UN-dtt))
        if(expt.lt.0.d0) then
**        write(6,*) ' problemme dans viscoe expt negatif ',expt
          expt=0.d0
        endif
        dpot=((EPSE-ED0-(UN/At)*expt**(UN/Bt))/XMdt)**XNdt
        dpot=MAX( dpot, XZERO)

        DTTDT=dtt+dpot*DELTAT
        DTTDT=MIN ( DTTDT , UN-1.D-05 )
**      dpoC=((EPSE-ED0-(UN/Ac)*(dct/(UN-dct))**(UN/Bc))/XMdC)**XNdc
        expc=(dct/(UN-dct))
        if(expc.lt.0.d0) then
**        write(6,*) ' problemme dans viscoe expc negatif ',expc
          expc=0.d0
        endif
        dpoC=((EPSE-ED0-(UN/Ac)*expc**(UN/Bc))/XMdC)**XNdc
        dpoC=MAX( dpoC, XZERO)
        DCTDT=dct+dpoC*DELTAT
        DCTDT=MIN ( DCTDT , UN-1.D-05 )
C
        DTDT = ALFAT*DTTDT + ALFAC*DCTDT


C
C           on borne la valeur de D a 0.99999
C
        DTDT=MIN ( DTDT , UN-1.D-05 )
C
C           on teste la croissance de D
C
        DTDT=MAX ( DTDT , DT )
C
      RETURN
      END