SUBROUTINE UMAT (STRESS, STATEV, DDSDDE, SSE, SPD, SCD,
     1 RPL, DDSDDT, DRPLDE, DRPLDT,
     2 STRAN, DSTRAN, TIME, DTIME, TEMP, DTEMP, PREDEF, DPRED, CMNAME,
     3 NDI, NSHR, NTENS, NSTATV, PROPS, NPROPS, COORDS, DROT, PNEWDT,
     4 CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER, KSPT, JSTEP, KINC)
    
        INCLUDE 'ABA_PARAM.INC'
    
       CHARACTER*80 CMNAME, CPNAME
       DIMENSION STRESS(NTENS), STATEV(NSTATV),
     1 DDSDDE(NTENS, NTENS), DDSDDT(NTENS), DRPLDE(NTENS),
     2 STRAN(NTENS), DSTRAN(NTENS), TIME(2), PREDEF(1), DPRED(1),
     3 PROPS(NPROPS), COORDS(3), DROT(3, 3), DFGRD0(3, 3), DFGRD1(3, 3),
     4 JSTEP(4)
    
       DIMENSION CFULL(NTENS, NTENS), DELTAT(NTENS)
      DOUBLE PRECISION KK, SIGMA3C, TAU1C, TAU2C, G1C, G2C, G3C, ETA
      DOUBLE PRECISION DMG_OLD, DMG, BETA, BI, DELTA, DELTA1, DELTA2
       DOUBLE PRECISION DELTA_SHEAR, DELTA0_3, DELTA_0SHEAR
       DOUBLE PRECISION DELTA0, DELTAF, RT_OLD, RT, DELTA3
        DOUBLE PRECISION CF6, CF7, CF5, CF4, KSH, KB, TAU
    
       PARAMETER(ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0, THREE = 3.D0)
    
       KK = 100
      SIGMA3C = 26       
      TAU1C = 26
      TAU2C = 26
      G1C = 15
       G2C = 15
       G3C = 15
       ETA = 0.1
    
      DMG_OLD = STATEV(1)
    
    ! DISPLACEMENT AT THE END OF THE INCREMENT
      DO I = 1, NTENS
        DELTAT(I) = STRAN(I) + DSTRAN(I)
      END DO
    
      DELTA3 = DELTAT(1)
       DELTA1 = DELTAT(2)
    
       
      IF (NTENS .EQ. TWO) THEN
        DELTA2 = ZERO
      ELSE
        DELTA2 = DELTAT(3)
      END IF
    
       DELTA_SHEAR = SQRT(DELTAT(1)**2 + DELTA2**2)
    
    ! SHEAR MODE PENALTY STIFFNESS
      KSH = G1C / G2C * (TAU2C / SIGMA3C)**2 * KK
    
    ! MIXED-MODE RATIO WITH SHEAR PENALTY STIFFNESS
      CF6 = KSH * DELTA_SHEAR**2 + KK * MAX(DELTAT(3), ZERO)**2
      CF7 = KSH**2 * DELTA_SHEAR**2 + KK**2 * MAX(DELTAT(3), ZERO)**2
      BI = (KSH * DELTA_SHEAR**2) / CF6
    
    ! Mixed-mode penalty stiffness
       KB = KK * (ONE - BI) + KSH * BI
    
    ! DISPLACEMENT JUMP NORM
      DELTA = CF6 / SQRT(CF7)
    
    ! PURE MODE ONSET DISPLACEMENT JUMP
       DELTA_0SHEAR = TAU1C / KSH
      DELTA0_3 = SIGMA3C / KK
    
    ! MIXED-MODE ONSET DISPLACEMENT JUMP
       CF4 = (KSH * DELTA_0SHEAR**2 - KK * DELTA0_3**2) * (BI**ETA)
        DELTA0 = SQRT((KK * DELTA0_3**2 + CF4) / KB)
    
    ! MIXED-MODE FINAL DISPLACEMENT JUMP
       DELTAC_SHEAR = TWO * G2C / TAU1C
       DELTAC_3 = TWO * G1C / SIGMA3C
      CF5 = (KSH * DELTA_0SHEAR * DELTAC_SHEAR - KK * DELTA0_3 *
     & DELTAC_3)
       DELTAF = (KK * DELTA0_3 * DELTAC_3 + CF5 * 
     & (BI**ETA)) / (KB * DELTA0)
    
    ! INTERNAL VARIABLE
      RT_OLD = (DELTA0 * DELTAF) / (DELTAF - DMG_OLD *
     & (DELTAF - DELTA0))
      RT = MAX(RT_OLD, DELTA)
    
      IF (DELTA .GT. RT_OLD) THEN
        DMG = (DELTAF * (RT - DELTA0)) / (RT * (DELTAF - DELTA0))
          DMG = MIN(DMG, ONE)
        ELSE
        DMG = DMG_OLD
      END IF
    
    ! FULL STIFFNESS MATRIX
      DO I = 1, NTENS
        DO J = 1, NTENS
            CFULL(I, J) = ZERO
        END DO
      END DO
    
       IF (DELTA3 .GE. ZERO) THEN
        CFULL(1, 1) = KK * (ONE - DMG)
       ELSE
        CFULL(1, 1) = KK
       END IF
    
       CFULL(2, 2) = KSH * (ONE - DMG)
    
    ! CFULL(3, 3) ONLY EXISTS IN 3D
      IF (NTENS .EQ. THREE) THEN
        CFULL(3, 3) = KSH * (ONE - DMG)
       END IF
    
    ! STRESS CALCULATION
      DO I = 1, NTENS
        STRESS(I) = ZERO
        DO J = 1, NTENS
            STRESS(I) = STRESS(I) + CFULL(I, J) * DELTAT(J)
        END DO
      END DO
    
    ! TANGENT STIFFNESS MATRIX
      DO I = 1, NTENS
        DO J = 1, NTENS
            DDSDDE(I, J) = ZERO
        END DO
      END DO
    
      IF (DELTA .GT. RT_OLD .AND. DELTA .LT. DELTAF) THEN
        COEFF = ((DELTAF * DELTA0) / ((DELTAF - DELTA0) * DELTA**3))
        
        IF (DELTA3 .GE. ZERO) THEN
            DDSDDE(1, 1) = CFULL(1, 1) - KK * COEFF * DELTA3**2
        ELSE
            DDSDDE(1, 1) = KK
        END IF
        
        DDSDDE(2, 2) = CFULL(2, 2) - KSH * COEFF * DELTA1**2
        DDSDDE(3, 3) = CFULL(3, 3) - KSH * COEFF * DELTA2**2
        DDSDDE(1, 2) = -KK * COEFF * DELTA1 * DELTA3
        DDSDDE(2, 1) = -KSH * COEFF * DELTA1 * DELTA3
        
        IF (NTENS .EQ. THREE) THEN
            DDSDDE(1, 3) = -KK * COEFF * DELTA2 * DELTA3
            DDSDDE(3, 1) = -KSH * COEFF * DELTA2 * DELTA3
            DDSDDE(2, 3) = -KSH * COEFF * DELTA1 * DELTA2
            DDSDDE(3, 2) = DDSDDE(2, 3)
        END IF
       ELSE
        DDSDDE(1, 1) = CFULL(1, 1)
        DDSDDE(2, 2) = CFULL(2, 2)
        
        IF (NTENS .EQ. THREE) THEN
            DDSDDE(3, 3) = CFULL(3, 3)
        END IF
       END IF
    
      STATEV(1) = DMG
      RETURN
      END