C 
      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 EELAS(6),EPLAS(6),DSTRES(6),PSTRES(6),DDFDDS(6),ATRT(6),
	1DDGDDS(6),DDGDS(6),DEVIAS(6),NBLFMT(6),R0(6),NBLQR(6),RB(6),IS(6),
     2JS(6),JCOBFZ(6,6),DDBDDS(6,6),DDEDDB(6,6),Q(6,6),R(6,6)
C      
	DIMENSION TEMP1(6,6),TEMP2(6,6),DDSBDDS(6,6),DDJ3DDS(6,6),
	1DSBDS(6),DJ3DS(6),DTHETDS(6),DAFADS(6),DC2MCDS(6),DC3MCDS(6),
	2DC2DS(6),DC3DS(6)
C
	REAL INVAR2,INVAR3,YIELDD,DDFDDJ,FB,FC,NBLFM1,NBLFM2,
	1DNBL,JCOBFM,NBLFZ,PHI,PUSI,NBLFMT,NBLFM,JCOBFZ,
     2DSTRES,PSTRES,DDFDDS,ATRT,DDGDDS,DDGDS,DEVIAS,RB,DDBDDS,
	3DDEDDB,R,Q,NBLQR,R0
C      
	REAL  PI,SIGMAM,TEMP1,TEMP2,
	1KF,DKF,KG,DKG,DDKG

	INTEGER NVALUE,IS,JS,L
C
      PARAMETER(ZERO=0.E0,ONE=1.E0,TWO=2.E0,THREE=3.E0,FOUR=4.E0,
	1SIX=6.E0,ENUMAX=0.499E0,NEWTON=20,TOLER=1.0E-6,PI=3.1415927)
C PROPS(1)-----E
C PROPS(2)-----NU
C PROPS(3)-----PHI
C PROPS(4)-----PUSI,DILATION ANGLE
C PROPS(5..)---SYIELD AND HARDING DATA
C CALLS AHARD FOR CURVE OF YIELD STRESS VS. PLASTIC STRAIN
C ************************************************
C ************************************************     
C 
	NBLFZ=ZERO
	NBLFM1=ZERO
	NBLFM2=ZERO
	JCOBFM=ZERO
C   
      EMOD=PROPS(1)
      ENU=PROPS(2)
	PHI=PROPS(3)
	PUSI=PROPS(4)
C 3K
      EBULK3=EMOD/(ONE-TWO*ENU)
C 2G
      EG2=EMOD/(ONE+ENU)
      EG=EG2/TWO
      EG3=THREE*EG
      ELAM=(EBULK3-EG2)/THREE
C   
      DO K1=1,NTENS
		DO K2=1,NTENS
			DDSDDE(K1,K2)=ZERO
			JCOBFZ(K1,K2)=ZERO
              DDBDDS(K1,K2)=ZERO
		END DO
		DSTRES(K1)=ZERO                   
	    DDGDS(K1)=ZERO
		ATRT(K1)=ZERO
		RB(K1)=ZERO
		NBLFMT(K1)=ZERO
      END DO
C ******************** 
      DO K1=1,NDI
		DO K2=1,NDI
			DDSDDE(K2,K1)=ELAM
	    END DO
          DDSDDE(K1,K1)=EG2+ELAM
	END DO
	DO K1=NDI+1,NTENS
		DDSDDE(K1,K1)=EG
	END DO
C 
	CALL ROTSIG(STATEV(      1),DROT,EELAS,2,NDI,NSHR)
	CALL ROTSIG(STATEV(NTENS+1),DROT,EPLAS,2,NDI,NSHR)
C 
	EQPLAS=STATEV(1+2*NTENS)
C **************************************************************************
C **************************************************************************
C **************************************************************************
C **************************************************************************
	DO K1=1,NTENS
		DO K2=1,NTENS
			DSTRES(K2)=DSTRES(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
		END DO
		EELAS(K1)=EELAS(K1)+DSTRAN(K1)
	END DO
	DO K1=1,NTENS
		PSTRES(K1)=STRESS(K1)+DSTRES(K1)
	END DO
C  
	SIGMAM=(PSTRES(1)+PSTRES(2)+PSTRES(3))/THREE
	DO K1=1,NDI
	 DEVIAS(K1)=PSTRES(K1)-SIGMAM
	END DO	
	INVAR2=(DEVIAS(1)**2+DEVIAS(2)**2+DEVIAS(3)**2)/TWO
	DO K1=NDI+1,NTENS
		INVAR2=INVAR2+PSTRES(K1)**2
	END DO
      SBAR=SQRT(INVAR2)
      INVAR3=DEVIAS(1)*DEVIAS(2)*DEVIAS(3)+
	1      2*PSTRES(4)*PSTRES(5)*PSTRES(6)-DEVIAS(1)*PSTRES(6)**2-
     2      DEVIAS(2)*PSTRES(5)**2-DEVIAS(3)*PSTRES(4)**2
	S3TA=-1.5*SQRT(THREE)*INVAR3/(SBAR**3)
      IF (S3TA.GT.1.0) THEN
	  S3TA=1.0
	ELSE IF (S3TA.LT.-1.0) THEN
        S3TA=-1.0
	END IF
      THETA=ASIN(S3TA)/THREE
      THETAT=(25.0/180.0)*PI
      A1=1.432052062044E0
      A2=0.406941858374E0
      B1=0.544290524902E0
	B2=0.673903324498E0
	IF (ABS(THETA).LE.THETAT) THEN
	 KF=COS(THETA)-SIN(PHI)*SIN(THETA)/SQRT(THREE)
	ELSE
       SGN=SIGN(ONE,THETA)
       CONST1=A1+A2*SGN*SIN(PHI)
	 CONST2=B1*SGN+B2*SIN(PHI)
       KF=CONST1-CONST2*S3TA
	END IF
C    
	NVALUE=(NPROPS-4)/2
	CALL AHARD(SYIEL0,HARD,EQPLAS,PROPS(5),NVALUE) 
      YIELDD=SYIEL0
      PARAMA=0.05*YIELDD
	FB=SIGMAM*SIN(PHI)+
	1   SQRT(SBAR*SBAR*KF*KF+PARAMA*PARAMA*COS(PHI)*COS(PHI))-
     2   YIELDD*COS(PHI)
C *********************************************************************
C *********************************************************************
C *********************************************************************
C *********************************************************************
	IF(FB.GT.TOLER) THEN
C*****
         ALFA=SBAR*KF
         ALFA=ALFA/SQRT(ALFA*ALFA+PARAMA*PARAMA*COS(PHI)*COS(PHI))	   
	  IF (ABS(THETA).LE.THETAT) THEN
	   DKF=-SIN(THETA)-SIN(PHI)*COS(THETA)/SQRT(THREE)
	  ELSE
         DKF=-THREE*CONST2*COS(THREE*THETA)
	  END IF
	   C1=SIN(PHI)
	   C2=KF-TAN(THREE*THETA)*DKF
	   C2=ALFA*C2
 	  IF (ABS(THETA).LE.THETAT) THEN
 	   C2=KF-TAN(THREE*THETA)*DKF
	   C2=ALFA*C2
         C3=-SQRT(THREE)*DKF/(TWO*COS(THREE*THETA)*SBAR*SBAR)
	   C3=ALFA*C3 
	  ELSE
 	   C2=KF
	   C2=ALFA*C2
         C3=ZERO
	  END IF
C	
	  DDFDDS(1)=C1/THREE+C2*DEVIAS(1)/(2*SBAR)+
	1    C3*(DEVIAS(2)*DEVIAS(3)-PSTRES(6)*PSTRES(6)+SBAR*SBAR/THREE)
	  DDFDDS(2)=C1/THREE+C2*DEVIAS(2)/(2*SBAR)+
	1    C3*(DEVIAS(1)*DEVIAS(3)-PSTRES(5)*PSTRES(5)+SBAR*SBAR/THREE)	  
	  DDFDDS(3)=C1/THREE+C2*DEVIAS(3)/(2*SBAR)+
	1    C3*(DEVIAS(1)*DEVIAS(2)-PSTRES(4)*PSTRES(4)+SBAR*SBAR/THREE)	  
	  DDFDDS(4)=C2*PSTRES(4)/SBAR+
     1           C3*TWO*(PSTRES(6)*PSTRES(5)-DEVIAS(3)*PSTRES(4))	  
	  DDFDDS(5)=C2*PSTRES(5)/SBAR+
     1           C3*TWO*(PSTRES(4)*PSTRES(6)-DEVIAS(2)*PSTRES(5))
	  DDFDDS(6)=C2*PSTRES(6)/SBAR+
     1           C3*TWO*(PSTRES(5)*PSTRES(4)-DEVIAS(1)*PSTRES(6))
C*****
	  PARAMB=0.2*YIELDD
	  IF (ABS(THETA).LT.THETAT) THEN
         KG=COS(THETA)-SIN(PUSI)*SIN(THETA)/SQRT(THREE)
	  ELSE
         SGN=SIGN(ONE,THETA)
         CONSTA=A1+A2*SGN*SIN(PUSI)
	   CONSTB=B1*SGN+B2*SIN(PUSI)
         KG=CONSTA-CONSTB*S3TA	
	  END IF
	   ALFAG=SBAR*KG
         ALFAG=ALFAG/SQRT(ALFAG*ALFAG+PARAMB*PARAMB*COS(PUSI)*COS(PUSI))	   
	  IF (ABS(THETA).LE.THETAT) THEN
	   DKG=-SIN(THETA)-SIN(PUSI)*COS(THETA)/SQRT(THREE)
	  ELSE
         DKG=-THREE*CONSTB*COS(THREE*THETA)
	  END IF
	   GC1=SIN(PUSI)
 	  IF (ABS(THETA).LE.THETAT) THEN
	   GC2=KG-TAN(THREE*THETA)*DKG
	   GC2=ALFAG*GC2
         GC3=-SQRT(THREE)*DKG/(TWO*COS(THREE*THETA)*SBAR*SBAR)
	   GC3=ALFAG*GC3 
	  ELSE
	   GC2=KG
	   GC2=ALFAG*GC2
	   GC3=ZERO 
	  END IF
C	
	  DDGDDS(1)=GC1/THREE+GC2*DEVIAS(1)/(2*SBAR)+
	1    GC3*(DEVIAS(2)*DEVIAS(3)-PSTRES(6)*PSTRES(6)+SBAR*SBAR/THREE)
	  DDGDDS(2)=GC1/THREE+GC2*DEVIAS(2)/(2*SBAR)+
	1    GC3*(DEVIAS(1)*DEVIAS(3)-PSTRES(5)*PSTRES(5)+SBAR*SBAR/THREE)	  
	  DDGDDS(3)=GC1/THREE+GC2*DEVIAS(3)/(2*SBAR)+
	1    GC3*(DEVIAS(1)*DEVIAS(2)-PSTRES(4)*PSTRES(4)+SBAR*SBAR/THREE)	  
	  DDGDDS(4)=GC2*PSTRES(4)/SBAR+
     1           GC3*TWO*(PSTRES(6)*PSTRES(5)-DEVIAS(3)*PSTRES(4))	  
	  DDGDDS(5)=GC2*PSTRES(5)/SBAR+
     1           GC3*TWO*(PSTRES(4)*PSTRES(6)-DEVIAS(2)*PSTRES(5))
	  DDGDDS(6)=GC2*PSTRES(6)/SBAR+
     1           GC3*TWO*(PSTRES(5)*PSTRES(4)-DEVIAS(1)*PSTRES(6))
C*****
			DO K1=1,NTENS
				DO K2=1,NTENS
					DDGDS(K2)=DDGDS(K2)+DDSDDE(K2,K1)*DDGDDS(K1)
				END DO
			END DO
C******
			B=(DDGDDS(1)-DDGDDS(2))**2+
	1			(DDGDDS(2)-DDGDDS(3))**2+(DDGDDS(3)-DDGDDS(1))**2
			DO K1=NDI+1,NTENS
				B=B+SIX*DDGDDS(K1)**2
			END DO
			B=SQRT(TWO*B)/THREE
			NBLFM2=HARD*B
C*****
			DO K1=1,NTENS
				NBLFM1=NBLFM1+DDFDDS(K1)*DDGDS(K1)
			END DO	
C******
			DNBL=FB/(NBLFM1+NBLFM2)
C**********************************************************************
C**********************************************************************
			DO K1=1,NTENS
				STRESS(K1)=PSTRES(K1)-DNBL*DDGDS(K1)
			END DO
C
C**********************************************************************
C**********************************************************************
	 DO KNEWTON=1,NEWTON
C***
	  SIGMAM=(STRESS(1)+STRESS(2)+STRESS(3))/THREE
	  DO K1=1,NDI
	   DEVIAS(K1)=STRESS(K1)-SIGMAM
	  END DO	
	  INVAR2=(DEVIAS(1)**2+DEVIAS(2)**2+DEVIAS(3)**2)/TWO
	  DO K1=NDI+1,NTENS
		INVAR2=INVAR2+STRESS(K1)**2
	  END DO
        SBAR=SQRT(INVAR2)
        INVAR3=DEVIAS(1)*DEVIAS(2)*DEVIAS(3)+
	1      2*STRESS(4)*STRESS(5)*STRESS(6)-DEVIAS(1)*STRESS(6)**2-
     2      DEVIAS(2)*STRESS(5)**2-DEVIAS(3)*STRESS(4)**2
	  S3TA=-1.5*SQRT(THREE)*INVAR3/(SBAR**3)
        IF (S3TA.GT.1.0) THEN
	   S3TA=1.0
	  ELSE IF (S3TA.LT.-1.0) THEN
         S3TA=-1.0
	  END IF
        THETA=ASIN(S3TA)/THREE
        THETAT=(25.0/180.0)*PI
        A1=1.432052062044E0
        A2=0.406941858374E0
        B1=0.544290524902E0
	  B2=0.673903324498E0
	  IF (ABS(THETA).LE.THETAT) THEN
	   KF=COS(THETA)-SIN(PHI)*SIN(THETA)/SQRT(THREE)
	  ELSE
         SGN=SIGN(ONE,THETA)
         CONST1=A1+A2*SGN*SIN(PHI)
	   CONST2=B1*SGN+B2*SIN(PHI)
         KF=CONST1-CONST2*S3TA
	  END IF
C
        NVALUE=(NPROPS-4)/2
 	  CALL AHARD(SYIEL0,HARD,EQPLAS,PROPS(5),NVALUE) 
        YIELDD=SYIEL0
	  PARAMA=0.05*YIELDD
C*****
         ALFA=SBAR*KF
         ALFA=ALFA/SQRT(ALFA*ALFA+PARAMA*PARAMA*COS(PHI)*COS(PHI))	   
	  IF (ABS(THETA).LE.THETAT) THEN
	   DKF=-SIN(THETA)-SIN(PHI)*COS(THETA)/SQRT(THREE)
	  ELSE
         DKF=-THREE*CONST2*COS(THREE*THETA)
	  END IF
	   C1=SIN(PHI)
	  IF (ABS(THETA).LE.THETAT) THEN
	   C2=KF-TAN(THREE*THETA)*DKF
	   C2=ALFA*C2
         C3=-SQRT(THREE)*DKF/(TWO*COS(THREE*THETA)*SBAR*SBAR)
	   C3=ALFA*C3 
	  ELSE
	   C2=KF
	   C2=ALFA*C2
	   C3=ZERO
	  END IF
C	
	  DDFDDS(1)=C1/THREE+C2*DEVIAS(1)/(2*SBAR)+
	1    C3*(DEVIAS(2)*DEVIAS(3)-STRESS(6)*STRESS(6)+SBAR*SBAR/THREE)
	  DDFDDS(2)=C1/THREE+C2*DEVIAS(2)/(2*SBAR)+
	1    C3*(DEVIAS(1)*DEVIAS(3)-STRESS(5)*STRESS(5)+SBAR*SBAR/THREE)	  
	  DDFDDS(3)=C1/THREE+C2*DEVIAS(3)/(2*SBAR)+
	1    C3*(DEVIAS(1)*DEVIAS(2)-STRESS(4)*STRESS(4)+SBAR*SBAR/THREE)	  
	  DDFDDS(4)=C2*STRESS(4)/SBAR+
     1           C3*TWO*(STRESS(6)*STRESS(5)-DEVIAS(3)*STRESS(4))	  
	  DDFDDS(5)=C2*STRESS(5)/SBAR+
     1           C3*TWO*(STRESS(4)*STRESS(6)-DEVIAS(2)*STRESS(5))
	  DDFDDS(6)=C2*STRESS(6)/SBAR+
     1           C3*TWO*(STRESS(5)*STRESS(4)-DEVIAS(1)*STRESS(6))
C******
  	  PARAMB=0.2*YIELDD
	  IF (ABS(THETA).LE.THETAT) THEN
         KG=COS(THETA)-SIN(PUSI)*SIN(THETA)/SQRT(THREE)
	  ELSE
         SGN=SIGN(ONE,THETA)
         CONSTA=A1+A2*SGN*SIN(PUSI)
	   CONSTB=B1*SGN+B2*SIN(PUSI)
         KG=CONSTA-CONSTB*S3TA	
	  END IF
	   ALFAG=SBAR*KG
         ALFAG=ALFAG/SQRT(ALFAG*ALFAG+PARAMB*PARAMB*COS(PUSI)*COS(PUSI))	   
	  IF (ABS(THETA).LE.THETAT) THEN
	   DKG=-SIN(THETA)-SIN(PUSI)*COS(THETA)/SQRT(THREE)
	  ELSE
         DKG=-THREE*CONSTB*COS(THREE*THETA)
	  END IF
	   GC1=SIN(PUSI)
	  IF (ABS(THETA).LE.THETAT) THEN
 	   GC2=KG-TAN(THREE*THETA)*DKG
	   GC2=ALFAG*GC2
         GC3=-SQRT(THREE)*DKG/(TWO*COS(THREE*THETA)*SBAR*SBAR)
	   GC3=ALFAG*GC3 
	  ELSE
	   GC2=KG
	   GC2=ALFAG*GC2
	   GC3=ZERO
	  END IF
C	
	  DDGDDS(1)=GC1/THREE+GC2*DEVIAS(1)/(2*SBAR)+
	1    GC3*(DEVIAS(2)*DEVIAS(3)-STRESS(6)*STRESS(6)+SBAR*SBAR/THREE)
	  DDGDDS(2)=GC1/THREE+GC2*DEVIAS(2)/(2*SBAR)+
	1    GC3*(DEVIAS(1)*DEVIAS(3)-STRESS(5)*STRESS(5)+SBAR*SBAR/THREE)	  
	  DDGDDS(3)=GC1/THREE+GC2*DEVIAS(3)/(2*SBAR)+
	1    GC3*(DEVIAS(1)*DEVIAS(2)-STRESS(4)*STRESS(4)+SBAR*SBAR/THREE)	  
	  DDGDDS(4)=GC2*STRESS(4)/SBAR+
     1           GC3*TWO*(STRESS(6)*STRESS(5)-DEVIAS(3)*STRESS(4))	  
	  DDGDDS(5)=GC2*STRESS(5)/SBAR+
     1           GC3*TWO*(STRESS(4)*STRESS(6)-DEVIAS(2)*STRESS(5))
	  DDGDDS(6)=GC2*STRESS(6)/SBAR+
     1           GC3*TWO*(STRESS(5)*STRESS(4)-DEVIAS(1)*STRESS(6))
C***
			DO K1=1,NTENS
				DDGDS(K1)=ZERO
			END DO
			DO K1=1,NTENS
				DO K2=1,NTENS
					DDGDS(K2)=DDGDS(K2)+DDSDDE(K2,K1)*DDGDDS(K1)
				END DO
			END DO
C*****
			B=(DDGDDS(1)-DDGDDS(2))**2+
	1			(DDGDDS(2)-DDGDDS(3))**2+(DDGDDS(3)-DDGDDS(1))**2
			DO K1=NDI+1,NTENS
				B=B+SIX*DDGDDS(K1)**2
			END DO
			B=SQRT(TWO*B)/THREE
			EQPLAS=STATEV(13)+DNBL*B
			CALL AHARD(SYIEL0,HARD,EQPLAS,PROPS(5),NVALUE)
			NBLFM2=HARD*B
C	
		 YIELDD=SYIEL0
           PARAMA=0.05*YIELDD/TAN(PHI)
 	     FC=SIGMAM*SIN(PHI)+
	1        SQRT(SBAR*SBAR*KF*KF+PARAMA*PARAMA*COS(PHI)*COS(PHI))-
     2        YIELDD*COS(PHI)					
C**************************************************************
	   IF (ABS(THETA).LE.THETAT) THEN
	    DDKG=-COS(THETA)+SIN(PUSI)*SIN(THETA)/SQRT(THREE)
	   ELSE
          DDKG=9.0*CONSTB*SIN(THREE*THETA)
	   END IF
         DO K1=1,NDI
	    DSBDS(K1)=DEVIAS(K1)/(TWO*SBAR)
         END DO
	   DO K1=NDI+1,NTENS
	    DSBDS(K1)=STRESS(K1)/SBAR
	   END DO
	  DJ3DS(1)=DEVIAS(2)*DEVIAS(3)-STRESS(6)*STRESS(6)+SBAR*SBAR/THREE
	  DJ3DS(2)=DEVIAS(1)*DEVIAS(3)-STRESS(5)*STRESS(5)+SBAR*SBAR/THREE
	  DJ3DS(3)=DEVIAS(1)*DEVIAS(2)-STRESS(4)*STRESS(4)+SBAR*SBAR/THREE
	  DJ3DS(4)=TWO*(STRESS(6)*STRESS(5)-DEVIAS(3)*STRESS(4))
	  DJ3DS(5)=TWO*(STRESS(4)*STRESS(6)-DEVIAS(2)*STRESS(5))
	  DJ3DS(6)=TWO*(STRESS(5)*STRESS(4)-DEVIAS(1)*STRESS(6))
C
	  IF (ABS(THETA).LE.THETAT) THEN
         DO K1=1,NTENS
         DTHETDS(K1)=SQRT(THREE)*(DJ3DS(K1)-THREE*INVAR3*DSBDS(K1)/SBAR)
         DTHETDS(K1)=-DTHETDS(K1)/(TWO*SBAR*SBAR*SBAR*COS(THREE*THETA))
	   END DO
	  ELSE
         DO K1=1,NTENS	  	
          DTHETDS(K1)=ZERO	
	   END DO	
	  END IF	
C	
	  DO K1=1,NTENS
         TERM1=ONE-ALFAG*ALFAG
         TERM2=SQRT(SBAR*SBAR*KG*KG+PARAMB*PARAMB*COS(PUSI)*COS(PUSI))
	   DAFADS(K1)=TERM1*(DSBDS(K1)*KG+SBAR*DKG*DTHETDS(K1))/TERM2
	  END DO
C
	  IF (ABS(THETA).LE.THETAT) THEN
         C2MC=KG-TAN(THREE*THETA)*DKG
         C3MC=-SQRT(THREE)*DKG/(TWO*COS(THREE*THETA)*SBAR*SBAR)
	  ELSE
         C2MC=KG
         C3MC=ZERO
	  END IF

	  IF (ABS(THETA).LE.THETAT) THEN
         DO K1=1,NTENS
          DC2MCDS(K1)=DKG+DDKG*TAN(THREE*THETA)-
	1              THREE*DKG/(COS(THREE*THETA)*COS(THREE*THETA))
          DC2MCDS(K1)=DC2MCDS(K1)*DTHETDS(K1)	
	   END DO
	  ELSE
          DC2MCDS(K1)=ZERO
	  END IF
C
	  IF (ABS(THETA).LE.THETAT) THEN
         DO K1=1,NTENS
          DC3MCDS(K1)=(DDKG-THREE*DKG*TAN(THREE*THETA))*DTHETDS(K1)+
	1              TWO*DKG*DSBDS(K1)/SBAR       
          DC3MCDS(K1)=DC3MCDS(K1)*SQRT(THREE)	
          DC3MCDS(K1)=DC3MCDS(K1)/(TWO*SBAR*SBAR*SBAR*COS(THREE*THETA))	
	   END DO
	  ELSE
          DC3MCDS(K1)=ZERO
	  END IF
C
        DO K1=1,NTENS
         DC2DS(K1)=ALFAG*DC2MCDS(K1)+C2MC*DAFADS(K1)
         DC3DS(K1)=ALFAG*DC3MCDS(K1)+C3MC*DAFADS(K1)
	  END DO
C		 
	  DO K1=1,NTENS
         TEMP1(1,K1)=DC2DS(1)*DSBDS(K1)
         TEMP1(2,K1)=DC2DS(2)*DSBDS(K1)
         TEMP1(3,K1)=DC2DS(3)*DSBDS(K1)
         TEMP1(4,K1)=DC2DS(4)*DSBDS(K1)
         TEMP1(5,K1)=DC2DS(5)*DSBDS(K1)
         TEMP1(6,K1)=DC2DS(6)*DSBDS(K1)		 
	  END DO
	  DO K1=1,NTENS
         TEMP2(1,K1)=DC3DS(1)*DJ3DS(K1)
         TEMP2(2,K1)=DC3DS(2)*DJ3DS(K1)
         TEMP2(3,K1)=DC3DS(3)*DJ3DS(K1)
         TEMP2(4,K1)=DC3DS(4)*DJ3DS(K1)
         TEMP2(5,K1)=DC3DS(5)*DJ3DS(K1)
         TEMP2(6,K1)=DC3DS(6)*DJ3DS(K1)		 
	  END DO
C
	  DO K1=1,NDI
         DDSBDDS(K1,K1)=ONE/THREE-DEVIAS(K1)*DEVIAS(K1)/(FOUR*SBAR*SBAR)
	  END DO
	  DO K1=NDI+1,NTENS
         DDSBDDS(K1,K1)=ONE-STRESS(K1)*STRESS(K1)/(SBAR*SBAR)
	  END DO
         DDSBDDS(2,1)=-ONE/SIX-DEVIAS(1)*DEVIAS(2)/(FOUR*SBAR*SBAR)
         DDSBDDS(3,1)=-ONE/SIX-DEVIAS(1)*DEVIAS(3)/(FOUR*SBAR*SBAR)
         DDSBDDS(3,2)=-ONE/SIX-DEVIAS(2)*DEVIAS(3)/(FOUR*SBAR*SBAR)
         DDSBDDS(4,1)=-STRESS(4)*DEVIAS(1)/(TWO*SBAR*SBAR)
         DDSBDDS(4,2)=-STRESS(4)*DEVIAS(2)/(TWO*SBAR*SBAR)
         DDSBDDS(4,3)=-STRESS(4)*DEVIAS(3)/(TWO*SBAR*SBAR)
         DDSBDDS(5,1)=-STRESS(5)*DEVIAS(1)/(TWO*SBAR*SBAR)
         DDSBDDS(5,2)=-STRESS(5)*DEVIAS(2)/(TWO*SBAR*SBAR)
         DDSBDDS(5,3)=-STRESS(5)*DEVIAS(3)/(TWO*SBAR*SBAR)
         DDSBDDS(5,4)=-STRESS(5)*DEVIAS(4)/(TWO*SBAR*SBAR)
         DDSBDDS(6,1)=-STRESS(6)*DEVIAS(1)/(TWO*SBAR*SBAR)
         DDSBDDS(6,2)=-STRESS(6)*DEVIAS(2)/(TWO*SBAR*SBAR)
         DDSBDDS(6,3)=-STRESS(6)*DEVIAS(3)/(TWO*SBAR*SBAR)
         DDSBDDS(6,4)=-STRESS(6)*DEVIAS(4)/(TWO*SBAR*SBAR)
         DDSBDDS(6,5)=-STRESS(6)*DEVIAS(5)/(TWO*SBAR*SBAR)
         DDSBDDS(1,2)=DDSBDDS(2,1)
         DDSBDDS(1,3)=DDSBDDS(3,1)
         DDSBDDS(2,3)=DDSBDDS(3,2)
         DDSBDDS(1,4)=DDSBDDS(4,1)
         DDSBDDS(2,4)=DDSBDDS(4,2)
         DDSBDDS(3,4)=DDSBDDS(4,3)
         DDSBDDS(1,5)=DDSBDDS(5,1)
         DDSBDDS(2,5)=DDSBDDS(5,2)
         DDSBDDS(3,5)=DDSBDDS(5,3)
         DDSBDDS(4,5)=DDSBDDS(5,4)
         DDSBDDS(1,6)=DDSBDDS(6,1)
         DDSBDDS(2,6)=DDSBDDS(6,2)
         DDSBDDS(3,6)=DDSBDDS(6,3)
         DDSBDDS(4,6)=DDSBDDS(6,4)
         DDSBDDS(5,6)=DDSBDDS(6,5)
C
         DDJ3DDS(1,1)=DEVIAS(1)-DEVIAS(2)-DEVIAS(3)
         DDJ3DDS(2,2)=DEVIAS(2)-DEVIAS(1)-DEVIAS(3)
         DDJ3DDS(3,3)=DEVIAS(3)-DEVIAS(1)-DEVIAS(2)
         DDJ3DDS(4,4)=-SIX*DEVIAS(3)
         DDJ3DDS(5,5)=-SIX*DEVIAS(2)
         DDJ3DDS(6,6)=-SIX*DEVIAS(1)
         DDJ3DDS(2,1)=TWO*DEVIAS(3)
         DDJ3DDS(3,1)=TWO*DEVIAS(2)
         DDJ3DDS(3,2)=TWO*DEVIAS(1)
         DDJ3DDS(4,1)=TWO*STRESS(4)
         DDJ3DDS(4,2)=TWO*STRESS(4)
         DDJ3DDS(4,3)=-FOUR*STRESS(4)
         DDJ3DDS(5,1)=TWO*STRESS(5)
         DDJ3DDS(5,2)=-FOUR*STRESS(5)
         DDJ3DDS(5,3)=TWO*STRESS(5)
         DDJ3DDS(5,4)=SIX*STRESS(6)
         DDJ3DDS(6,1)=-FOUR*STRESS(6)
         DDJ3DDS(6,2)=TWO*STRESS(6)
         DDJ3DDS(6,3)=TWO*STRESS(6)
         DDJ3DDS(6,4)=SIX*STRESS(5)
         DDJ3DDS(6,5)=SIX*STRESS(4)        
	   DDJ3DDS(1,2)=DDJ3DDS(2,1)
         DDJ3DDS(1,3)=DDJ3DDS(3,1)
         DDJ3DDS(2,3)=DDJ3DDS(3,2)
         DDJ3DDS(1,4)=DDJ3DDS(4,1)
         DDJ3DDS(2,4)=DDJ3DDS(4,2)
         DDJ3DDS(3,4)=DDJ3DDS(4,3)
         DDJ3DDS(1,5)=DDJ3DDS(5,1)
         DDJ3DDS(2,5)=DDJ3DDS(5,2)
         DDJ3DDS(3,5)=DDJ3DDS(5,3)
         DDJ3DDS(4,5)=DDJ3DDS(5,4)
         DDJ3DDS(1,6)=DDJ3DDS(6,1)
         DDJ3DDS(2,6)=DDJ3DDS(6,2)
         DDJ3DDS(3,6)=DDJ3DDS(6,3)
         DDJ3DDS(4,6)=DDJ3DDS(6,4)
         DDJ3DDS(5,6)=DDJ3DDS(6,5)
C
	   DO K1=1,NTENS
		  DO K2=1,NTENS
			DDSBDDS(K1,K2)=DDSBDDS(K1,K2)/SBAR
			DDJ3DDS(K1,K2)=DDJ3DDS(K1,K2)/THREE
		  END DO
	   END DO
C
	   DO K1=1,NTENS
		  DO K2=1,NTENS
              DDBDDS(K1,K2)=TEMP1(K1,K2)+GC2*DDSBDDS(K1,K2)+
	1                 TEMP2(K1,K2)+GC3*DDJ3DDS(K1,K2)   
		  END DO
	   END DO
C       
			DO 50 K1=1,NTENS
			DO 50 K2=1,NTENS
				DDEDDB(K1,K2)=ZERO
				DO 10 K3=1,NTENS
					DDEDDB(K1,K2)=DDEDDB(K1,K2)+
	1					DDSDDE(K1,K3)*DDBDDS(K3,K2)
10				CONTINUE
50			CONTINUE
C******
			DO K1=1,NTENS
				DO K2=1,NTENS
					Q(K1,K2)=DNBL*DDEDDB(K1,K2)
				END DO
			END DO
			DO K1=1,NTENS
				Q(K1,K1)=Q(K1,K1)+ONE
			END DO
C****
		CALL BRINV(Q,6,L,IS,JS)
C
			IF(ABS(FC).LT.10.0) GO TO 100
C
C**************************************************************
			DO K1=1,NTENS
				NBLFMT(K1)=ZERO
			END DO
			DO K1=1,NTENS
				DO K2=1,NTENS
					NBLFMT(K2)=NBLFMT(K2)+Q(K2,K1)*DDGDS(K1)	                   			
				END DO
			END DO
			NBLFM1=ZERO	
			DO K1=1,NTENS
				NBLFM1=NBLFM1+NBLFMT(K1)*DDFDDS(K1)
			END DO
			NBLFM=NBLFM1+NBLFM2
C******
			DO K1=1,NTENS
				R0(K1)=STRESS(K1)-PSTRES(K1)+DNBL*DDGDS(K1)
			END DO
C*****
			DO K1=1,NTENS
				NBLQR(K1)=ZERO
			END DO
			DO K1=1,NTENS
				DO K2=1,NTENS
					NBLQR(K2)=NBLQR(K2)+Q(K2,K1)*R0(K1)	
				END DO
			END DO
			NBLFZ=ZERO
			DO K1=1,NTENS
				NBLFZ=NBLFZ+DDFDDS(K1)*NBLQR(K1)
			END DO
			NBLFZ=FC-NBLFZ
C*****
			CNBL=NBLFZ/NBLFM
C*****
			DO K1=1,NTENS
				DSTRES(K1)=-NBLQR(K1)-CNBL*NBLFMT(K1)
			END DO
			DO K1=1,NTENS
				STRESS(K1)=STRESS(K1)+DSTRES(K1)
			END DO	
C****
			DNBL=DNBL+CNBL
	 END DO
C
C***
		WRITE(7,2) NEWTON,FC
2			FORMAT(//,'***WARNING-PLASTICITY LAGORITHM DID NOT',
	1		'CONVERGE AFTER',I3,' ITERRATIONS',2X,'FC EQUALS',E11.4)
C****
100		CONTINUE
C
C*****
		DO 53 K1=1,NTENS
			EPLAS(K1)=EPLAS(K1)+DNBL*DDGDDS(K1)
			EELAS(K1)=EELAS(K1)-DNBL*DDGDDS(K1)
53		CONTINUE
C
C******** 
			DO 150 K1=1,NTENS
			DO 150 K2=1,NTENS
				R(K1,K2)=ZERO
				DO 110 K3=1,NTENS
					R(K1,K2)=R(K1,K2)+Q(K1,K3)*DDSDDE(K3,K2)
110				CONTINUE
150			CONTINUE
C*****
			DO K1=1,NTENS
				DO K2=1,NTENS
					ATRT(K2)=ATRT(K2)+R(K2,K1)*DDFDDS(K1)
					RB(K2)=RB(K2)+R(K2,K1)*DDGDDS(K1)
				END DO
			END DO
C
			DO K1=1,NTENS
				DO K2=1,NTENS
					JCOBFZ(K1,K2)=RB(K1)*ATRT(K2)
				END DO
			END DO
C*****
			DO K1=1,NTENS
				JCOBFM=JCOBFM+DDFDDS(K1)/RB(K1)
			END DO
			JCOBFM=JCOBFM-NBLFM2
C*****
			DO K1=1,NTENS
				DO K2=1,NTENS
					DDSDDE(K1,K2)=R(K1,K2)-JCOBFZ(K1,K2)/JCOBFM
				END DO 
			END DO
C************
	ELSE
		DO K1=1,NTENS
			STRESS(K1)=PSTRES(K1)
		END DO	
C
	END IF
C***********
	DO K1=1,NTENS
		STATEV(K1)=EELAS(K1)		
		STATEV(K1+NTENS)=EPLAS(K1)	
	END DO
	STATEV(1+2*NTENS)=EQPLAS
C
C
C
C
	RETURN
	END
C
*****************************************************************************
      SUBROUTINE AHARD(SYIELD,HARD,EQPLAS,TABLE,NVALUE)
C
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION TABLE(2,NVALUE)
C
C    SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO
      SYIELD=TABLE(1,NVALUE)
      HARD=0.0
C
C   IF MORE THAN ONE ENTRY, SEARCH TABLE
C
      IF(NVALUE.GT.1) THEN
        DO 10 K1=1,NVALUE-1
           EQPL1=TABLE(2,K1+1)
           IF(EQPLAS.LT.EQPL1) THEN
             EQPL0=TABLE(2,K1)
             IF(EQPL1.LE.EQPL0) THEN
                WRITE(6,1)
 1              FORMAT(//,30X,'***ERROR - PLASTIC STRAIN MUST BE ',
     1                 'ENTERED IN ASCENDING ORDER')
                CALL XIT
              ENDIF
C
C           CURRENT YIELD STRESS AND HARDENING
C
            DEQPL=EQPL1-EQPL0
            SYIEL0=TABLE(1,K1)
            SYIEL1=TABLE(1,K1+1)
            DSYIEL=SYIEL1-SYIEL0
            HARD=DSYIEL/DEQPL
            SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD
            GOTO 20
            ENDIF
 10         CONTINUE
 20         CONTINUE
      ENDIF
      RETURN
      END
C
C******
	SUBROUTINE BRINV(Q,N,L,IS,JS)
	DIMENSION Q(N,N),IS(N),JS(N)
	REAL Q,T,D
	L=1
	DO 300 K=1,N
		D=0.0
		DO 310 I=K,N
		DO 310 J=K,N
			IF(ABS(Q(I,J)).GT.D) THEN
				D=ABS(Q(I,J))
				IS(K)=I
				JS(K)=J
			END IF
310		CONTINUE
		IF(D+1.0.EQ.1.0)THEN
			L=0
C			WRITE(*,320)
			RETURN 
		END IF
C320		FORMAT(1X,'ERR**NOT INV')
		DO 330 J=1,N
			T=Q(K,J)
			Q(K,J)=Q(IS(K),J)
			Q(IS(K),J)=T
330		CONTINUE
		DO 340 I=1,N
			T=Q(I,K)
			Q(I,K)=Q(I,JS(K))
			Q(I,JS(K))=T
340		CONTINUE
		Q(K,K)=1/Q(K,K)
		DO 350 J=1,N
			IF(J.NE.K)THEN
				Q(K,J)=Q(K,J)*Q(K,K)
			END IF
350		CONTINUE
		DO 370 I=1,N
			IF(I.NE.K)THEN
				DO 360 J=1,N
					IF(J.NE.K)THEN
						Q(I,J)=Q(I,J)-Q(I,K)*Q(K,J)
					END IF
360				CONTINUE
			END IF
370		CONTINUE
		DO 380 I=1,N
			IF(I.NE.K)THEN
				Q(I,K)=-Q(I,K)*Q(K,K)
			END IF
380		CONTINUE
300	CONTINUE
	DO 430 K=N,1,-1
		DO 410 J=1,N
			T=Q(K,J)
			Q(K,J)=Q(JS(K),J)
			Q(JS(K),J)=T
410		CONTINUE
		DO 420 I=1,N
			T=Q(I,K)
			Q(I,K)=Q(I,IS(K))
			Q(I,IS(K))=T
420		CONTINUE
430	CONTINUE
	RETURN
	END
C
C