       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,KSTEP,KINC)
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*8 CMNAME
      DIMENSION STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),
     1 DDSDDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS),
     2 TIME(2),PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3),
     3 DFGRD0(3,3),DFGRD1(3,3),EPLAS(6),DEPLAS(6),ALPHA(6),DALPHA(6),
     4 DEVS(6),FLOW(6),OLDS(6),DSTRESS(6),ALPHA1(6),ALPHA2(6),ALPHA3(6),
     5 DALPHA1(6), DALPHA2(6), DALPHA3(6), ZETA(6), FLOWEPS(6),
     6 DZETA(6),DEVOPE(NTENS,NTENS),SUMH(NTENS,NTENS),H1(NTENS,NTENS),
     7 H2(NTENS,NTENS),H3(NTENS,NTENS),ZL(NTENS,NTENS),
     8 ZINVL(NTENS,NTENS),ZB(NTENS,NTENS+NTENS)
      PARAMETER (ZERO=0.0, ONE=1.0, NEWTONMAX=10000,
     1 TWO=2.0, THREE=3.0, TOLER=0.0001, TEN=10.0, FOUR=4.0)
C
      EMOD=PROPS(1)
      ENU=PROPS(2)
      YIELD=PROPS(3)
      C1=PROPS(4)
      GAMMA1=PROPS(5)
      C2=PROPS(6)
      GAMMA2=PROPS(7)
      C3=PROPS(8)
      GAMMA3=PROPS(9)
      AQ=PROPS(10)
      BQ=PROPS(11)
      CQ=PROPS(12)
      B=PROPS(13)
      H=PROPS(14)
      ETA=PROPS(15)
C
      EBULK3=EMOD/(ONE-TWO*ENU)
      EG=EMOD/((ONE+ENU)*TWO)
      ELAM=(EBULK3-EG*TWO)/THREE
C
C ELASTIC STIFFNESS
      DO K1=1, NDI
       DO K2=1, NDI
        DDSDDE(K2, K1)=ELAM
       END DO
       DDSDDE(K1, K1)=EG*TWO+ELAM
      END DO
      DO K1=NDI+1, NTENS
       DDSDDE(K1 ,K1)=EG
      END DO
C
      DO K1=1, NTENS
       EPLAS(K1)=STATEV(K1)
       ALPHA1(K1)=STATEV(K1+NTENS)
       ALPHA2(K1)=STATEV(K1+2*NTENS)
       ALPHA3(K1)=STATEV(K1+3*NTENS)
       ALPHA(K1)=ALPHA1(K1)+ALPHA2(K1)+ALPHA3(K1)
       ZETA(K1)=STATEV(K1+4*NTENS)
      END DO
      QMEMORY=STATEV(5*NTENS+1)
      R1=STATEV(5*NTENS+2)
      R2=STATEV(5*NTENS+3)
      P=STATEV(5*NTENS+4)
C
C ELASTICKY PREDIKTOR
      DO K1=1, NTENS
       OLDS(K1)=STRESS(K1)
      END DO
      DO K1=1, NTENS
       DO K2=1, NTENS
        STRESS(K1)=STRESS(K1)+DDSDDE(K1, K2)*DSTRAN(K2)
       END DO
      END DO
C
C MISSESOVO EKV NAPETI
      SMISESTR=(STRESS(1)-ALPHA(1)-STRESS(2)+ALPHA(2))**TWO
     1 +(STRESS(2)-ALPHA(2)-STRESS(3)+ALPHA(3))**TWO
     2 +(STRESS(3)-ALPHA(3)-STRESS(1)+ALPHA(1))**TWO
      DO K1=NDI+1,NTENS
       SMISESTR=SMISESTR+6.*(STRESS(K1)-ALPHA(K1))**TWO
      END DO
      SMISESTR=SQRT(SMISESTR/TWO)
C
C DEVIATOR NAPETI PREDIKTOR
      SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE
      DO K1=1, NDI
       DEVS(K1)=STRESS(K1)-SHYDRO
      END DO
      DO K1=NDI+1, NTENS
       DEVS(K1)=STRESS(K1)
      END DO
C
C NORMALA
      DO K1=1,NTENS
       FLOW(K1)=THREE/TWO*(DEVS(K1)-ALPHA(K1))/SMISESTR
      END DO
C
C PLASTICITY FUNCTION
      F=SMISESTR-YIELD-R1-R2
      DO K1=1, NTENS
       DSTRESS(K1)=STRESS(K1)-OLDS(K1)
      END DO
C PROMITNUTI PRIRUSTKU DO NORMALY
      PRODUCT=ZERO
C
      DO K1=1, NDI
       PRODUCT=PRODUCT+DSTRESS(K1)*FLOW(K1)
      END DO
      DO K1=NDI+1, NTENS
       PRODUCT=PRODUCT+TWO*DSTRESS(K1)*FLOW(K1)
      END DO
      DO K1=1, NTENS
       DALPHA(K1)=ZERO
       DALPHA1(K1)=ZERO
       DALPHA2(K1)=ZERO
       DALPHA3(K1)=ZERO
       DEPLAS(K1)=ZERO
       DZETA(K1)=ZERO
      END DO
C
      DEPL=ZERO
      DDEPL=ZERO
      DEPL_OLD=ZERO
      DR1=ZERO
      DR2=ZERO
      DQMEMORY=ZERO
      Q=AQ+BQ*(ONE-EXP(-CQ*QMEMORY))
C
      IF(F.LE.ZERO) GOTO 20
      IF(PRODUCT.LE.ZERO) GOTO 20
C
C ITERACE
      DO NEWTON=1, NEWTONMAX
C Y
       Y=B*(Q-R1)/((ONE+B*DEPL)**TWO)+H*(ONE+B*P+TWO*B*DEPL)
C DR
       DR1=B*(Q-R1)*DEPL/(ONE+B*DEPL)
       DR2=H*(ONE+B*(P+DEPL))*DEPL
C Z
       Z=ZERO
       DO K1=1, NDI
        Z=Z+(DEVS(K1)-ALPHA(K1)-DALPHA(K1))*(((TWO/THREE*C1*FLOW(K1)
     1 -GAMMA1*ALPHA1(K1))*(ONE+GAMMA1*DEPL)-(TWO/THREE*C1*DEPL*FLOW(K1)
     2 -GAMMA1*ALPHA1(K1)*DEPL)*GAMMA1)/((ONE+GAMMA1*DEPL)**TWO)+((TWO
     3 /THREE*C2*FLOW(K1)-GAMMA2*ALPHA2(K1))*(ONE+GAMMA2*DEPL)-(TWO
     4 /THREE*C2*DEPL*FLOW(K1)-GAMMA2*ALPHA2(K1)*DEPL)*GAMMA2)/((ONE
     5 +GAMMA2*DEPL)**TWO)+((TWO/THREE*C3*FLOW(K1)-GAMMA3*ALPHA3(K1))
     6 *(ONE+GAMMA3*DEPL)-(TWO/THREE*C3*DEPL*FLOW(K1)-GAMMA3*ALPHA3(K1)
     7 *DEPL)*GAMMA3)/((ONE+GAMMA3*DEPL)**TWO))
       END DO
       DO K1=NDI+1, NTENS
        Z=Z+TWO*(DEVS(K1)-ALPHA(K1)-DALPHA(K1))*(((TWO/THREE*C1*FLOW(K1)
     1 -GAMMA1*ALPHA1(K1))*(ONE+GAMMA1*DEPL)-(TWO/THREE*C1*DEPL*FLOW(K1)
     2 -GAMMA1*ALPHA1(K1)*DEPL)*GAMMA1)/((ONE+GAMMA1*DEPL)**TWO)+((TWO
     3 /THREE*C2*FLOW(K1)-GAMMA2*ALPHA2(K1))*(ONE+GAMMA2*DEPL)-(TWO
     4 /THREE*C2*DEPL*FLOW(K1)-GAMMA2*ALPHA2(K1)*DEPL)*GAMMA2)/((ONE
     5 +GAMMA2*DEPL)**TWO)+((TWO/THREE*C3*FLOW(K1)-GAMMA3*ALPHA3(K1))
     5 *(ONE+GAMMA3*DEPL)-(TWO/THREE*C3*DEPL*FLOW(K1)-GAMMA3*ALPHA3(K1)
     7 *DEPL)*GAMMA3)/((ONE+GAMMA3*DEPL)**TWO))
       END DO
       Z=Z*((-THREE/TWO)/SMISESTR)
C DDEPL A DEPL
       DDEPL=(SMISESTR-THREE*EG*DEPL-YIELD-R1-R2-DR1-DR2)/(THREE*EG-Z+Y)
       DEPL=DEPL+DDEPL
C P MA BYT KLADNE
       IF (DEPL .LT. ZERO) THEN
        DEPL=DEPL_OLD/TEN
       END IF
C DALPHA
       DO K1=1, NTENS
        DALPHA1(K1)=(TWO/THREE*C1*DEPL*FLOW(K1)-GAMMA1*ALPHA1(K1)
     1 *DEPL)/(ONE+GAMMA1*DEPL)
        DALPHA2(K1)=(TWO/THREE*C2*DEPL*FLOW(K1)-GAMMA2*ALPHA2(K1)
     1 *DEPL)/(ONE+GAMMA2*DEPL)
        DALPHA3(K1)=(TWO/THREE*C3*DEPL*FLOW(K1)-GAMMA3*ALPHA3(K1)
     1 *DEPL)/(ONE+GAMMA3*DEPL)
        DALPHA(K1)=DALPHA1(K1)+DALPHA2(K1)+DALPHA3(K1)
       END DO
C SMISESTR
       SMISESTR=(STRESS(1)-ALPHA(1)-DALPHA(1)-STRESS(2)+ALPHA(2)
     1 +DALPHA(2))**TWO+(STRESS(2)-ALPHA(2)-DALPHA(2)-STRESS(3)
     2 +ALPHA(3)+DALPHA(3))**TWO+(STRESS(3)-ALPHA(3)-DALPHA(3)
     3 -STRESS(1)+ALPHA(1)+DALPHA(1))**TWO
       DO K1=NDI+1,NTENS
        SMISESTR=SMISESTR+6.*(STRESS(K1)-ALPHA(K1)-DALPHA(K1))**TWO
       END DO
       SMISESTR=SQRT(SMISESTR/TWO)
C FLOW
       DO K1=1, NTENS
        FLOW(K1)=THREE/TWO*(DEVS(K1)-ALPHA(K1)-DALPHA(K1))/SMISESTR
       END DO
       XKRITERIUM=ABS(ONE-DEPL_OLD/DEPL)
       DEPL_OLD=DEPL
       IF (XKRITERIUM .LE. TOLER) GOTO 10
       IF(NEWTON.GT.(NEWTONMAX-TWO)) CALL XIT
      END DO
C
10    CONTINUE
C PLASTICKA DEFORMACE A PLASTICKY KOREKTOR
      DO K1=1, NDI
       DEPLAS(K1)=DEPL*FLOW(K1)
      END DO
      DO K1=NDI+1, NTENS
       DEPLAS(K1)=TWO*DEPL*FLOW(K1)
      END DO
      DO K1=1, NTENS
       STRESS(K1)=OLDS(K1)
       DO K2=1, NTENS
        STRESS(K1)=STRESS(K1)+DDSDDE(K1, K2)*(DSTRAN(K2)-DEPLAS(K2))
       END DO
      END DO
C
C UPDATE PAMETOVE PLOCHY
C EKVIVALENTNI PLASTICKA DEFORMACE
      SMISESE=ZERO
      DO K1=1,NDI
       SMISESE=SMISESE+(EPLAS(K1)-ZETA(K1))**TWO
      END DO
      DO K1=NDI+1,NTENS
       SMISESE=SMISESE+TWO*(EPLAS(K1)/TWO-ZETA(K1))**TWO
      END DO
      SMISESE=SQRT(TWO/THREE*SMISESE)
      F=SMISESE-QMEMORY
      IF (F .LE. ZERO) GOTO 20
C NORMALA
      DO K1=1, NDI
       FLOWEPS(K1)=(EPLAS(K1)-ZETA(K1))/SMISESE
      END DO
      DO K1=NDI+1, NTENS
       FLOWEPS(K1)=(EPLAS(K1)/TWO-ZETA(K1))/SMISESE
      END DO
C ZMENY PAMETOVE PLOCHY
      DQMEMORY=ETA*F
      DO K1=1, NTENS
       DZETA(K1)=(ONE-ETA)*F*FLOWEPS(K1)
      END DO
C
C KONZISTETNI TECNY MODUL
      SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE
      DO K1=1, NDI
       DEVS(K1)=STRESS(K1)-SHYDRO
      END DO
      DO K1=NDI+1, NTENS
       DEVS(K1)=STRESS(K1)
      END DO
      DO K1=1, NTENS
       FLOW(K1)=THREE/TWO*(DEVS(K1)-ALPHA(K1)-DALPHA(K1))
     1 /(YIELD+R1+DR1+R2+DR2)
      END DO
C MATICE ID
      DO K1=1, NTENS
       DO K2=1, NTENS
       DEVOPE(K2,K1)=0
       END DO
      END DO
      DO K1=1, NDI
       DO K2=1, NDI
       DEVOPE(K2,K1)=-ONE/THREE
       END DO
       DEVOPE(K1,K1)=TWO/THREE
      END DO
      DO K1=NDI+1, NTENS
       DEVOPE(K1,K1)=ONE
      END DO
C MATICE  H
      DO K1=1, NTENS
       DO K2=1, NTENS
        H1(K1,K2)=-TWO/THREE*GAMMA1*(ALPHA1(K1)+DALPHA1(K1))*FLOW(K2)
        H2(K1,K2)=-TWO/THREE*GAMMA2*(ALPHA2(K1)+DALPHA2(K1))*FLOW(K2)
        H3(K1,K2)=-TWO/THREE*GAMMA3*(ALPHA3(K1)+DALPHA3(K1))*FLOW(K2)
       END DO
      END DO
      DO K1=1, NDI
        H1(K1,K1)=H1(K1,K1)+(TWO/THREE)*C1
        H2(K1,K1)=H2(K1,K1)+(TWO/THREE)*C2
        H3(K1,K1)=H3(K1,K1)+(TWO/THREE)*C3
      END DO
      DO K1=NDI+1, NTENS
        H1(K1,K1)=H1(K1,K1)+C1/THREE
        H2(K1,K1)=H2(K1,K1)+C2/THREE
        H3(K1,K1)=H3(K1,K1)+C3/THREE
      END DO
      DO K1=1, NTENS
       DO K2=1, NTENS
        H1(K1,K2)=H1(K1,K2)/(ONE+GAMMA1*DEPL)
        H2(K1,K2)=H2(K1,K2)/(ONE+GAMMA2*DEPL)
        H3(K1,K2)=H3(K1,K2)/(ONE+GAMMA3*DEPL)
        SUMH(K1,K2)=H1(K1,K2)+H2(K1,K2)+H3(K1,K2)
       END DO
      END DO
C MATICE L
      TEMP=B*(Q-R1-DR1)/(ONE+B*DEPL)+H*(ONE+B*P+TWO*B*DEPL)
      DO K1=1, NDI
       DO K2=1, NTENS
        ZL(K1,K2)=SUMH(K1,K2)+(TWO/THREE)*FLOW(K1)*FLOW(K2)
     1 *(TEMP-(YIELD+R1+R2+DR1+DR2)/DEPL)*(TWO/THREE)
       END DO
      ZL(K1,K1)=ZL(K1,K1)+TWO/THREE*(YIELD+R1+R2+DR1+DR2)/DEPL+TWO*EG
      END DO
      DO K1=NDI+1, NTENS
       DO K2=1, NTENS
        ZL(K1,K2)=TWO*SUMH(K1,K2)+TWO*(TWO/THREE)*FLOW(K1)*FLOW(K2)
     1 *(TEMP-(YIELD+R1+R2+DR1+DR2)/DEPL)*(TWO/THREE)
       END DO
      ZL(K1,K1)=ZL(K1,K1)+TWO/THREE*(YIELD+R1+R2+DR1+DR2)/DEPL+TWO*EG
      END DO
C INVERZE L
C
C Download URL: http://wp.me/p61TQ-zb
C Last modified: 2010/12/30
C A general purpose matrix inverter by augmenting-pivoting technique:
C A B C | 1 0 0 1 0 0 | J K L
C D E F | 0 1 0 => 0 1 0 | M N O
C G H I | 0 0 1 0 0 1 | P Q R
C Based on a lecture by Prof. McFarland
C http://math.uww.edu/~mcfarlat/inverse.htm
C Explanation of passed parameters:
C
C MAKE AUGMENTED MATRIX
      DO I=1,NTENS
       DO J=1,NTENS
        ZB(I,J)=0.0
        ZB(I,J+NTENS)=0.0
        ZB(I,J)=ZL(I,J)
        IF(I.EQ.J) THEN
         ZB(I,J+NTENS)=1.0
        END IF
       END DO
      END DO
C
      DO I=1,NTENS
C CHOOSE THE LEFTMOST NON-ZERO ELEMENT AS PIVOT
       DO J=1,NTENS
        IF(ABS(ZB(I,J)).GT.0.0)THEN
         PIVOT=ZB(I,J)
         EXIT
        END IF
       END DO
C STEP 1: Change the chosen pivot into "1" by dividing
C the pivot's row by the pivot number
       DO J=1,2*NTENS
        ZB(I,J)=ZB(I,J)/PIVOT
       END DO
       PIVOT=ZB(I,I) !UPDATE PIVOT VALUE
C STEP 2: Change the remainder of the pivot's COLUMN into 0's
C by adding to each row a suitable multiple of the PIVOT ROW
       DO K=1,NTENS !ROW
        IF(K.NE.I) THEN
         XNUM=ZB(K,I)/PIVOT !SAME COLUMN WITH THE CURRENT PIVOT
         DO J=1,2*NTENS !COL
          ZB(K,J)=ZB(K,J)-XNUM*ZB(I,J)
         END DO
        END IF
       END DO
      END DO
C PREPARE THE FINAL INVERTED MATRIX
      DO I=1,NTENS
       DO J=1,NTENS
        ZINVL(I,J)=ZB(I,J+NTENS)
       END DO
      END DO
C
C L KRAT M2
      DO K1=NDI+1, NTENS
       DO K2=1, NTENS
        ZINVL(K1,K2)=ZINVL(K1,K2)/TWO
       END DO
      END DO
C TANGENCIALNI MATICE
      DO K1=1, NTENS
       DO K2=1, NTENS
        DO K3=1, NTENS
        TEMP=FOUR*EG*EG*ZINVL(K1,K3)*DEVOPE(K3,K2)
        DDSDDE(K1,K2)=DDSDDE(K1,K2)-TEMP
        END DO
       END DO
      END DO
C
20    CONTINUE
C
C ULOZENI EPLAS, ALPHA A P
      DO K1=1, NTENS
       STATEV(K1)=EPLAS(K1)+DEPLAS(K1)
       STATEV(NTENS+K1)=ALPHA1(K1)+DALPHA1(K1)
       STATEV(2*NTENS+K1)=ALPHA2(K1)+DALPHA2(K1)
       STATEV(3*NTENS+K1)=ALPHA3(K1)+DALPHA3(K1)
       STATEV(4*NTENS+K1)=ZETA(K1)+DZETA(K1)
      END DO
      STATEV(5*NTENS+1)=QMEMORY+DQMEMORY
      STATEV(5*NTENS+2)=R1+DR1
      STATEV(5*NTENS+3)=R2+DR2
      STATEV(5*NTENS+4)=P+DEPL
C
      RETURN
      END
