C NMAZAR    SOURCE    LJASON    07/10/24    21:15:03     5930

       SUBROUTINE NMAZAR(XMAT,NMATT,VAR0,VARF,NVARI,SIG0,
     & SIGF,DEPST,NSTRS,ISTEP)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

* Entrees
       INTEGER NMATT, NVARI, NSTRS, ISTEP
       REAL*8 XMAT(NMATT), VAR0(NVARI), SIG0(NSTRS), DEPST(NSTRS)
       REAL*8 SIGF(NSTRS), VARF(NVARI)

* Variables intermediaires
       REAL*8 EPS33(3,3),EPSIPP(3),EPSILT(3),VALP33(3,3)
       REAL*8 SIGP(3),SIGPT(3),SIGPC(3),TRSIGT,TRSIGC
       REAL*8 SIG00(6), EPSTIM
       REAL*8 YOUN,XNU,EPSD0,ACOM,BCOM,ATRA,BTRA,BETA
       REAL*8 XZERO,UN,DEUX,XPETIT, DDHOOK(6,6)
       REAL*8 DINI,D,DT,DC,EPSTIL,ALFAT,ALFAC,GAMMA
       REAL*8 EPSILO(6), PREC
       INTEGER ISTRS,JSTRS,KCAS,IRTD
       CHARACTER CMATE
       PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12)

      CMATE = 'ISOTROPE'
       YOUN = XMAT(1)
       XNU  = XMAT(2)
       EPSD0= XMAT(9)
       ACOM = XMAT(5)
       BCOM = XMAT(7)
       ATRA = XMAT(6)
       BTRA = XMAT(8)
       BETA = XMAT(20)
       DINI = VAR0(2)

C
C     calcul des seuils
C
C
C        calcul de la deformation élastique
C
* Calcul de la matrice elastique
       KCAS=2
       CALL DOHMAS(XMAT,'ISOTROPE',2,NSTRS,KCAS,DDHOOK,IRTD)

       PREC=1.D-20
       CALL INVALM(DDHOOK,NSTRS,NSTRS,KERRE,PREC)

* Calcul de La partie élastique des déformations
       CALL MATVE1 (DDHOOK,SIG0,NSTRS,NSTRS,EPSILO,1)

*       DO ISTRS = 1,NSTRS
*         EPSILO(ISTRS) = 1/(1-DINI)*EPSILO(ISTRS) + DEPST(ISTRS)
*       END DO
C
C        calcul des deformations principales
C
C
C        on reecrit les deformations sous forme matricielle
C
      CALL ENDOCB (EPSILO,EPS33,2,2)
C
C        et on diagonalise
C
      CALL JACOB3 (EPS33,3,EPSIPP,VALP33)


      IF (ISTEP .EQ. 0 .OR. ISTEP.EQ.2) THEN
C
C       on calcule la matrice de hooke et les contraintes ppales
C
         CMATE = 'ISOTROPE'
         KCAS=2
         CALL DOHMAS(XMAT,'ISOTROPE',2,NSTRS,KCAS,DDHOOK,IRTD)
         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
      END IF

C        on calcule le epsilontild
C
      EPSTIL=MAX( XZERO , EPSIPP(1) )**2 +
     1       MAX( XZERO , EPSIPP(2) )**2 +
     2       MAX( XZERO , EPSIPP(3) )**2
      EPSTIL=SQRT (EPSTIL)

      IF (ISTEP .EQ. 0) THEN
         EPSTIM=EPSTIL
         VARF(1)=EPSTIL
      ELSE IF (ISTEP .EQ. 1) THEN
         VARF(9)=VAR0(9)
         VARF(3)=VAR0(3)
         VARF(4)=VAR0(4)
         VARF(5)=VAR0(5)
         VARF(6)=VAR0(6)
         VARF(7)=VAR0(7)
         VARF(8)=VAR0(8)
         VARF(2)=DINI
         VARF(1)=EPSTIL
      ELSE IF (ISTEP .EQ. 2) THEN
         EPSTIM=VAR0(1)
         VARF(1)=EPSTIM
      ELSE
         PRINT*,'DANS MAZARS ISTEP = 0,1,2 ET PAS ',ISTEP
      END IF

      IF ( (ISTEP .EQ. 0) .OR. (ISTEP .EQ. 2)) THEN
C
C      on calcule l'endommagement et les contraintes
C
C
C        on calcule ALFAT ALFAC DT et DC puis D
C        dans le cas ou le seuil initial est depasse
C
         IF ( EPSTIM .GT. EPSD0) THEN
C
C      calcul de l'endommagement
C
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/(EPSTIL*EPSTIL)
            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)+
*     1                                 SIGPC(3)*SIGPC(3)
*               GAMMA=-SQRT(GAMMA)/TRSIGC
*               EPSTIM=EPSTIM*GAMMA
*            END IF
C
C        amelioration de la reponse en cisaillement pour beta > 1.
C
            IF (BETA .GT. UN) THEN
               IF ( ALFAT .GT. XPETIT ) THEN
                  ALFAT=ALFAT**BETA
               END IF
               IF ( ALFAC .GT. XPETIT ) THEN
                  ALFAC=ALFAC**BETA
               END IF
            END IF

C
C        on calcule DT et DC puis D
C        dans le cas ou le seuil initial est depasse
C
C        on est oblige de verifier car on a pu multiplier par gamma
C
            IF (EPSTIM .GT. EPSD0) THEN
               DT=UN - EPSD0*(UN-ATRA)/EPSTIM -
     1            ATRA*EXP(-BTRA*(EPSTIM-EPSD0))
               DC=UN - EPSD0*(UN-ACOM)/EPSTIM -
     1            ACOM*EXP(-BCOM*(EPSTIM-EPSD0))
            ELSE
               DT=XZERO
               DC=XZERO
            END IF
            D = ALFAT*DT + ALFAC*DC
C
C           on borne la valeur de D a 0.99999
C
            D=MIN ( D , UN-1.D-05 )
         ELSE
            D=XZERO
         END IF
C
C           on teste la croissance de D
C
         D=MAX ( D , DINI )
C
C           on le stocke dans les variables internes finales
C
         VARF(2)= D
C
C           on calcule les contraintes finales
C
         DO 500 ISTRS=1,NSTRS
            SIGF(ISTRS)=SIG0(ISTRS)*(UN-D)
500      CONTINUE

*        Et on stocke la différence entre contraintes réelles et contraintes plastiques

         DO ISTRS=1,NSTRS
           VARF(ISTRS+2) = SIG0(ISTRS)-SIGF(ISTRS)
         END DO

      END IF

      RETURN
      END