      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,
     2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,
     3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,
     4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 MATERL
      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),
     4 DFGRD0(3,3),DFGRD1(3,3)
C
	DIMENSION PS(3),DSTRESS(NTENS)
	PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)

      EK=PROPS(1)
	EN=PROPS(2)
	RF=PROPS(3)
	C=PROPS(4)
	FAI=PROPS(5)/180.0*3.1415926
	DFAI=PROPS(6)/180.0*3.1415926
	EKB=PROPS(7)
	EM=PROPS(8)
	EKUR=PROPS(9)
	PA=PROPS(10)
	S1S3O=STATEV(1)
	S3O=STATEV(2)
	SSS=STATEV(3)
	
	CALL GETPS(STRESS,PS,NTENS)
	FAI=FAI-DFAI*LOG10(S3O/PA)
	CALL GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O,EKB,EM
     1	,SSS,S1S3O)

	EBULK3=EMOD/(ONE-TWO*ENU)
      EG2=EMOD/(ONE+ENU)
      EG=EG2/TWO
      EG3=THREE*EG
      ELAM=(EBULK3-EG2)/THREE
	CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)
	DSTRESS=0.0
	CALL GETSTRESS(DDSDDE,DSTRESS,DSTRAN,NTENS)
	DO 701 I1=1,NTENS
	STRESS(I1)=STRESS(I1)+DSTRESS(I1)
701	CONTINUE
	
	CALL GETPS(STRESS,PS,NTENS)
	
	CALL GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O,EKB,EM
     1	,SSS,S1S3O)


	EBULK3=EMOD/(ONE-TWO*ENU)
      EG2=EMOD/(ONE+ENU)
      EG=EG2/TWO
      EG3=THREE*EG
      ELAM=(EBULK3-EG2)/THREE
	
	CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)


	IF(PS(3).GT.S3O)S3O=PS(3)
 	IF((PS(1)-PS(3)).GT.S1S3O)S1S3O=PS(1)-PS(3)
	IF(S.GT.SSS)SSS=S
	STATEV(1)=S1S3O
	STATEV(2)=S3O
	STATEV(3)=SSS

	END
	


	SUBROUTINE GETPS(STRESS,PS,NTENS)
	INCLUDE 'ABA_PARAM.INC'
	DIMENSION PS(3),STRESS(NTENS)
	CALL SPRINC(STRESS,PS,1,3,3)
	DO 310 I=1,2
	DO 320 J=I+1,3
	IF(PS(I).GT.PS(J))THEN
	PPS=PS(I)
	PS(I)=PS(J)
	PS(J)=PPS
	END IF
320	CONTINUE
310	CONTINUE
	DO 330 K1=1,3
	PS(K1)=-PS(K1)
330	CONTINUE
	RETURN
	END 

	SUBROUTINE GETEMOD(PS,EK,EN,RF,C,FAI,ENU,PA,EKUR,EMOD,S,S3O
     1	,EKB,EM,SSS,S1S3O)

	INCLUDE 'ABA_PARAM.INC'
	DIMENSION PS(3)
	
  	S=(1-SIN(FAI))*(PS(1)-PS(3))
	IF(PS(3).LT.0) THEN
	PSFEI=0.1
	ELSE
	PSFEI=PS(3)
	END IF
	S=S/(2*C*COS(FAI)+2*PSFEI*SIN(FAI))
	IF(S.GE.0.99) THEN
	S=0.99
	END IF
	EMOD=EK*PA*((S3O/PA)**EN)*((1-RF*S)**2)
	IF((S.LT.SSS).AND.((PS(1)-PS(3)).LT.S1S3O))THEN
	EMOD=EKUR*PA*((S3O/PA)**EN)
	END IF

      AA=EKB*PA*((S3O/PA)**EMM)
	ENU=0.5*(1.0-EMOD/3./AA)
	IF(ENU.GT.0.49)ENU=0.49
	IF(ENU.LT.0.01)ENU=0.01
	END 


	SUBROUTINE GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)
	INCLUDE 'ABA_PARAM.INC'
	DIMENSION DDSDDE(NTENS,NTENS)
	DO 20 K1=1,NTENS
        DO 10 K2=1,NTENS
           DDSDDE(K2,K1)=0.0
 10     CONTINUE
 20   CONTINUE
      DO 40 K1=1,NDI
        DO 30 K2=1,NDI
           DDSDDE(K2,K1)=ELAM
 30     CONTINUE
        DDSDDE(K1,K1)=EG2+ELAM
 40   CONTINUE
      DO 50 K1=NDI+1,NTENS
        DDSDDE(K1,K1)=EG
 50   CONTINUE
	RETURN
	END

	SUBROUTINE GETSTRESS(DDSDDE,STRESS,DSTRAN,NTENS)
	INCLUDE 'ABA_PARAM.INC'
	DIMENSION DDSDDE(NTENS,NTENS),STRESS(NTENS),DSTRAN(NTENS)
	DO 70 K1=1,NTENS
        DO 60 K2=1,NTENS
           STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
 60     CONTINUE
 70   CONTINUE
	RETURN
	END
