      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*80 CMNAME
c
      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),IS(4),JS(4)
c*********************
C     Local Arrays   *
C*********************
	DIMENSION DA(4),DB(4),STRESSC(4),
     1PS(3),DDSDDET(4,4),DDS(4,4),DCB(4),STRESSB(4),
	2DSTRESS(4),DV(4),DAC(4),
	3DC(4),DR(4),DQMATRIX(4,4),DAQ(4),DQR(4),DQCB(4)
C
	DOUBLE PRECISION  DALFA,DBEIDA,PJ,DFAIL,DAQR,DAQCB,DGLING,DMIU,
     1DFRIC,DDILA,DLAMDAB,DLAMDAC,DTERM,DRR
	INTEGER KEY
c
	PARAMETER (ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,FOUR=4.D0,
     1SIX=6.D0,NINE=9.D0,EIGHT=8.D0,HALF=5.D-1,TOLER=1.D-4,NNEWTON=100)
C	
C**************************
C	ELASTIC PROPERTIES
C**************************
	DGLING=PROPS(1)
	DMIU=PROPS(2)
C**************************
C	MATERIAL PROPERTIES
C**************************
	DFRIC=PROPS(3)/180.D0*3.1415926
	DDILA=PROPS(4)/180.D0*3.1415926
C	EFRICV=PROPS(5)/180.0*3.1415926
c	SIGY0=PROPS(6)
C	WRITE (*,*) GLING,EMIU,SIGY0
C
C	ZERO MATRICES
C
	DO K1=1,4
	DO K2=1,4
		DDS(K1,K2)=ZERO
	END DO
	END DO
C
	DO K1=1,4
	DO K2=1,4
		DDSDDET(K1,K2)=ZERO
	END DO
	END DO
C
C
C
	EP=STATEV(1)
	WRITE (7,*) "NOEL=",NOEL
	WRITE (7,*) "NPT=",NPT
C
	DALFA=TWO/SQRT(THREE)*SIN(DFRIC)/(THREE-SIN(DFRIC))
	DBEIDA=TWO/SQRT(THREE)*SIN(DDILA)/(THREE-SIN(DDILA))
c	EGAMA=TWO/SQRT(THREE)*SIN(EFRICV)/(THREE-SIN(EFRICV))
C	WRITE (7,*) "*****************************************"
C	WRITE (7,*) DALFA,DBEIDA
C	WRITE (7,*) "*****************************************"
C
C	SET UP ELASTICITY MATRIX
C
	DBULK3=TWO*DGLING*(ONE+DMIU)/(ONE-TWO*DMIU)
	DLAM=(DBULK3-TWO*DGLING)/THREE
	DO K1=1,3
	DO K2=1,3
		DDS(K2,K1)=DLAM
	END DO
	   DDS(K1,K1)=DGLING*TWO+DLAM
	END DO
	   DDS(4,4)=DGLING

C
C	
C
      WRITE (7,*) "STRESS=",STRESS
	WRITE (7,*) "DSTRAN=",DSTRAN
	CALL KMLT1(DDS,DSTRAN,DV,NTENS)
	WRITE (7,*) "DV=",DV
	DO K1=1,4
		STRESSB(K1)=STRESS(K1)+DV(K1)
	END DO
	WRITE (7,*) "STRESSB=",STRESSB
c
	CALL SPRINC(STRESSB,PS,1,3,1)
	WRITE (7,*) "PS=",PS
C
	CALL YIELDSURFACE(STRESSB,DALFA,PJ,DFAIL)
	WRITE (7,*) "DFAIL=",DFAIL
c
      DLAMDAB=ZERO
	KEY=2
C
      IF (DFAIL.GT.ZERO) THEN
	WRITE (7,*) "ACTIVATE PLASTICITY"
	KEY=1
C	WRITE (7,*) "EZY,PJ,EYE1=",EZY,PJ,EALFA,EYE1
C	WRITE (7,*) "EYIELD=",EYIELD
C
C	IF (EFAIL.GT.(ONE+TOLER)*EYIELD) THEN
C	
C	CALCULATE THE PLASTIC FLOW DIRECTION
	WRITE (*,*) "ACTIVATE PLASTICITY"
C

C	write (7,*) "xndir=",xndir
C
C     
C
	CALL BFLOW(STRESSB,DALFA,DBEIDA,DA,DB,DC,PJ)
	WRITE (7,*) "DC=",DC
	WRITE (7,*) "DA=",DA
	WRITE (7,*) "DB=",DB
C
	CALL KMLT1(DDS,DB,DV,NTENS)
	WRITE (7,*) "DV=",DV
	CALL DOTPROD(DA,DV,DLAMDAB,NTENS)
	WRITE (7,*) "DA=",DA
	WRITE (7,*) "DV=",DV
	WRITE (7,*) "TERM=",DLAMDAB
	DLAMDAB=DFAIL/DLAMDAB
C	WRITE (7,*) "DLAMDAB,EZY,TERM2=",DLAMDAB,EZY,TERM2
C
C
C
C
C	WRITE (7,*) "DV=",DV
	DO K1=1,4
		DV(K1)=DLAMDAB*DV(K1)
	END DO
C	WRITE (7,*) "DSTRESS=",DSTRESS
	DO K1=1,4
		STRESSC(K1)=STRESSB(K1)-DV(K1)
	END DO
C
	DO 99 KNEWTON=1,NNEWTON
	WRITE (7,*) "THE START OF THE INTERATION NNEWTON"
	WRITE (7,*) knewton
C	WRITE (7,*) "STRESSC=",STRESSC
C	WRITE (7,*) "STRESSB=",STRESSB
C	WRITE (7,*) "DLAMDAB=",DLAMDAB
	CALL YIELDSURFACE(STRESSC,DYIELD,PJ,DFAIL)
	WRITE (*,*) "PJ,DFAIL AT C=",PJ,DFAIL
	CALL BFLOW(STRESSC,DALFA,DBEIDA,DA,DB,DC,PJ)
	WRITE (7,*) "DC=",DC
	WRITE (7,*) "DA=",DA
	WRITE (7,*) "DB=",DB
	CALL KMLT1(DDS,DB,DV,NTENS)
	WRITE (7,*) "DV=",DV
	WRITE (7,*) "DLAMDAB=",DLAMDAB
	DO K1=1,4
		DR(K1)=STRESSC(K1)-(STRESSB(K1)-DLAMDAB*DV(K1))
	END DO
	DRR=ZERO
	DO K1=1,4
		DRR=DRR+DR(K1)*DR(K1)
	END DO
		DRR=SQRT(DRR)
c		WRITE (7,*) "DR=",DR
		WRITE (7,*) "DRR=",DRR
	IF (DRR.LT.0.0001) THEN 
c	.AND.(ABS(ER(4)).LT.0.000001))
	GOTO 200
	ELSE
C
			IF (KNEWTON.GT.20) THEN
			WRITE (*,*) "AFTER",kNEWTON, "INTRATIONS, STOPPED"
			CALL EXIT
			ELSE
			WRITE (*,2) KNEWTON
			WRITE (*,*) "DRR=",DRR
2	FORMAT (//,'***WARNING -PLASTICITY ALGORITHM DID NOT CONVERGE',
     1            'AFTER',I4,'--ITERATIONS***') 
			END IF
C
	CALL PRODQM(PJ,DC,DLAMDAB,DDS,DQMATRIX)
	CALL INVERSE(DQMATRIX,L,IS,JS)
	CALL KMLT2(DA,DQMATRIX,DAQ,NTENS)
	CALL DOTPROD(DAQ,DR,DAQR,NTENS)
	WRITE (7,*) "DV=",DV
	CALL DOTPROD(DAQ,DV,DAQCB,NTENS)
	WRITE (7,*) "DFAIL,DAQR,DAQCB=",DFAIL,DAQR,DAQCB
	DLAMDAC=(DFAIL-DAQR)/DAQCB
	WRITE (7,*) "DLAMDAC=",DLAMDAC
	CALL KMLT1(DQMATRIX,DR,DQR,NTENS)
	CALL KMLT1(DQMATRIX,DV,DQCB,NTENS)
	DO K1=1,4
	STRESSC(K1)=STRESSC(K1)-DQR(K1)-DLAMDAC*DQCB(K1)
	END DO
	DLAMDAB=DLAMDAB+DLAMDAC
	END IF
99	CONTINUE
	END IF
C
C
C
C
200	WRITE (7,*) "MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM"
	WRITE (7,*) "MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM"
C
	IF (KEY.EQ.1) THEN
C
	DO K1=1,NTENS
	STRESSB(K1)=STRESSC(K1)
	END DO
C
	WRITE (7,*) "DB,DA=",DB,DA
	CALL KMLT1(DDS,DB,DCB,NTENS)
	CALL KMLT2(DA,DDS,DAC,NTENS)
	WRITE (7,*) "DA,DCB=",DA,DCB
	CALL DOTPROD(DA,DCB,DTERM,NTENS)
	WRITE (7,*) "DA,DCB=",DA,DCB
	WRITE (7,*) "DFACTOR=",DTERM
	CALL DYADICPROD(DAC,DCB,DDSDDET,NTENS)
C
	DO K1=1,4
	DO K2=1,4
		DDSDDET(K2,K1)=DDSDDET(K2,K1)/DTERM
	END DO
	END DO
C
	END IF
C	WRITE (7,*) "DDSDDET=",DDSDDET
C
	DO K1=1,NTENS
		STRESS(K1)=STRESSB(K1)
	END DO
	WRITE (7,*) "STRESS=",STRESS
C
	DO K1=1,4
	DO K2=1,4
		DDSDDE(K2,K1)=DDS(K2,K1)-DDSDDET(K2,K1)
	END DO
	END DO
	WRITE (7,*) "DDSDDE=",DDSDDE

C
	RETURN
	END 
C
C
C
C
C
C
C
c
C------------------------------------------------------------
C
C	SUBROUTINE TO CALCULATE MATRIXQ
C
C------------------------------------------------------------
C
	SUBROUTINE YIELDSURFACE(STR,DALF,DPJJ,DFAI)
C
	INCLUDE 'ABA_PARAM.INC'
C
	DIMENSION STR(4)
C
	PARAMETER (ZERO=0.D0, ONE=1.D0,TWO=2.D0,THREE=3.D0,FOUR=4.D0,
     1EIGHT=8.D0,NINE=9.D0,SIX=6.D0)
	DOUBLE PRECISION EALF, PJJ,EFAI

	DPJJ=ONE/SIX*((STR(1)-STR(2))**2+(STR(2)-STR(3))**2
     1+(STR(3)-STR(1))**2)+STR(4)*STR(4)
C
	DPJJ=SQRT(DPJJ)
	DFAI=DALF*(STR(1)+STR(2)+STR(3))+DPJJ
	WRITE (7,*) "DFAIL00=",DFAI
C
	RETURN
	END 
C
C
C
C
C
	SUBROUTINE BFLOW(STR,EALF,EBEID,DA,DB,DC,PJJ)
C
	INCLUDE 'ABA_PARAM.INC'
C
	DIMENSION STR(4),DA(4),DB(4),DC(4)
C
	DOUBLE PRECISION EALF, EBEID, PJJ, EYE
C
	PARAMETER (ZERO=0.D0, ONE=1.D0,TWO=2.D0,THREE=3.D0,FOUR=4.D0,
     1EIGHT=8.D0,NINE=9.D0)
C
	EYE=ONE/THREE*(STR(1)+STR(2)+STR(3))
C
      DO K = 1,3
      DC(K) = (ONE/TWO)*(STR(K)-EYE)/PJJ
      END DO
	DC(4) = STR(4)/PJJ
C
	DO K1=1,3
	DA(K1)=DC(K1)+EALF
	END DO
	DA(4)=DC(4)
	DO K1=1,3
	DB(K1)=DC(K1)+EBEID
	END DO
	DB(4)=DC(4)
	RETURN
	END 
C
C------------------------------------------------------------
C
C	SUBROUTINE TO CALCULATE MATRIXQ
C
C------------------------------------------------------------
C
	SUBROUTINE PRODQM(PJJ,DC,DLAMDA,DDS,DQMATRIX)
C
	INCLUDE 'ABA_PARAM.INC'
C
	DIMENSION DCC(4,4),DAA(4,4),DQMATRIX(4,4),DQMATRIXT(4,4),
     1DC(4),DIDEN(4,4),DDS(4,4)
	DOUBLE PRECISION PJJ, DLAMDA
C
	PARAMETER (ZERO=0.D0, ONE=1.D0,TWO=2.D0,THREE=3.D0,FOUR=4.D0,
     1EIGHT=8.D0,NINE=9.D0,SIX=6.D0)
C
	CALL DYADICPROD(DC,DC,DCC,NTENS)
	WRITE (*,*) "DCC=",DCC
C
	DO 20 I=1,4
      DO 20 J=1,4
        IF(I .EQ. J) THEN
          DIDEN(I,J)=1.0D0
        ELSE
          DIDEN(I,J)=0.0D0
        END IF
  20  CONTINUE
	
	DO 10 K1=1,3
		DO 10 K2=1,3
	IF (K1.EQ.K2) THEN
	DAA(K1,K2)=ONE/THREE
	ELSE 
	DAA(K1,K2)=-ONE/SIX
	END IF
10	CONTINUE
	DAA(4,4)=ONE
	WRITE (*,*) "DAA=", DAA
C
C	WRITE (*,*) "PJ=",PJ
	DO K1=1,4
	DO K2=1,4
	DQMATRIXT(K1,K2)=(DAA(K1,K2)-DCC(K1,K2))/PJJ
	END DO
	END DO
C	WRITE (*,*) "EQMATRIXT=",EQMATRIXT
C	WRITE (*,*) "DDS=",DDS
	CALL KMLT3(DDS,DQMATRIXT,DQMATRIX,4)
C	WRITE (*,*) "EQMATRIX1=", EQMATRIX

C	WRITE (*,*) "DLAMDAB=",DLAMDAB
	DO K1=1,4
	DO K2=1,4
	DQMATRIX(K1,K2)=DIDEN(K1,K2)+DLAMDA*DQMATRIX(K1,K2)
	END DO
	END DO
	WRITE (*,*) "DQMATRIX2=", DQMATRIX
	RETURN 
	END
**
**
**
**
***************************************************
**         MULTIPLY 4X4 MATRIX WITH 4X4 MATRIX    *
***************************************************
*USER SUBROUTINE
      SUBROUTINE KMLT3(DM1,DM2,DM,NTENS)
C      
      INCLUDE 'ABA_PARAM.INC'
C
      PARAMETER (M=4,ZERO=0.0)
C
      DIMENSION DM1(M,M),DM2(M,M),DM(M,M)
C
      DO 10 I=1,NTENS
      DO 10 J=1,NTENS 
	DM(I,J)=ZERO
	DO 20 K=1,NTENS
      DM(I,J)=DM(I,J)+DM1(I,K)*DM2(K,J)
20    CONTINUE
10    CONTINUE
      RETURN
      END
C
C
C
***************************************************
**         MULTIPLY 4X4 MATRIX WITH 4X1 VECTOR    *
***************************************************
*USER SUBROUTINE
      SUBROUTINE KMLT1(DM1,DM2,DM,NTENS)
C      
      INCLUDE 'ABA_PARAM.INC'
C
      PARAMETER (M=4)
C
      DIMENSION DM1(M,M),DM2(M),DM(M)
C
      DO 10 I=1,4
      X=0.0
      DO 20 K=1,4
      Y=DM1(I,K)*DM2(K)
      X=X+Y
20    CONTINUE
      DM(I)=X
10    CONTINUE
      RETURN
      END
**
**
***************************************************
**         MULTIPLY 1x4 MATRIX WITH 4X4 VECTOR    *
***************************************************
      SUBROUTINE KMLT2(DM1,DM2,DM,NTENS)
C      
      INCLUDE 'ABA_PARAM.INC'
C
      PARAMETER (M=4)
C
      DIMENSION DM1(M),DM2(M,M),DM(M)
C
      DO 10 I=1,NTENS
      X=0.0
      DO 20 K=1,NTENS 
      Y=DM1(K)*DM2(K,I)
      X=X+Y
20    CONTINUE
      DM(I)=X
10    CONTINUE
      RETURN
      END
**
********************************************
**         DOT PRODUCT OF TWO VECTORS      *
********************************************
*USER SUBROUTINE
      SUBROUTINE DOTPROD(DM1,DM2,DM,NTENS)
C      
      INCLUDE 'ABA_PARAM.INC'
C
      DIMENSION DM1(NTENS),DM2(NTENS)
C
      Y=0.D0
      DO 20 K=1,NTENS
      X=DM1(K)*DM2(K)
      Y=X+Y
20	CONTINUE
      DM=Y
      RETURN
      END
**
      SUBROUTINE DYADICPROD(DM1,DM2,DM3,NTENS)
C      
      INCLUDE 'ABA_PARAM.INC'
C
      PARAMETER (M=4)
C
      DIMENSION DM1(4),DM2(4),DM3(4,4)
C
      DO I=1,4
       DO J=1,4
        DM3(I,J) = DM1(I)*DM2(J)
       END DO
      END DO
C  
      RETURN
      END
**
**
C
C
C********************************************
**      CALCULATE INVERSE OF A MATRIX      *
********************************************
	SUBROUTINE INVERSE(EAA,L,IS,JS)
C
	INCLUDE 'ABA_PARAM.INC'
C
	PARAMETER (ZERO=0.D0, ONE=1.D0,TWO=2.D0,THREE=3.D0,FOUR=4.D0,
     1EIGHT=8.D0,NINE=9.D0)
C
	DIMENSION EAA(4,4),IS(4),JS(4)
	DOUBLE PRECISION EAA
	L=1
	DO 100 K=1,4
	   EDD=ZERO
		DO 10 I=K,4
		DO 10 J=K,4
			IF (ABS(EAA(I,J)).GT.EDD) THEN
				EDD=ABS(EAA(I,J))
				IS(K)=I
				JS(K)=J
			END IF
10	CONTINUE
C
	IF (EDD+ONE.EQ.ONE) THEN
		L=0
		WRITE(*,20)
		RETURN
	END IF
20	FORMAT (1X,'ERROR*****NOT INV')
C
	DO 30 J=1,4
		ETT=EAA(K,J)
		EAA(K,J)=EAA(IS(K),J)
		EAA(IS(K),J)=ETT
30	CONTINUE
C
	DO 40 I=1,4
		ETT=EAA(I,K)
		EAA(I,K)=EAA(I,JS(K))
		EAA(I,JS(K))=ETT
40	CONTINUE
C
	EAA(K,K)=ONE/EAA(K,K)
	DO 50 J=1,4
		IF (J.NE.K) THEN
			EAA(K,J)=EAA(K,J)*EAA(K,K)
		END IF
50	CONTINUE
C
	DO 70 I=1,4
		IF (I.NE.K)THEN
			DO 60 J=1,4
				IF (J.NE.K) THEN
					EAA(I,J)=EAA(I,J)-EAA(I,K)*EAA(K,J)
				END IF
60	        CONTINUE
		END IF
70	CONTINUE
C
	DO 80 I=1,4
		IF (I.NE.K) THEN
			EAA(I,K)=-EAA(I,K)*EAA(K,K)
		END IF
80	CONTINUE
100	CONTINUE
C
	DO 130 K=4,1,-1
		DO 110 J=1,4
			ETT=EAA(K,J)
			EAA(K,J)=EAA(JS(K),J)
			EAA(JS(K),J)=ETT
110		CONTINUE
		DO 120 I=1,4
			ETT=EAA(I,K)
			EAA(I,K)=EAA(I,IS(K))
			EAA(I,IS(K))=ETT
120		CONTINUE
130	CONTINUE
	RETURN
	END