SUBROUTINE UMAT(STRESS, STATEV, DDSDDE, SSE, SPC, SCD, RPL,DDSDDT,
     *                DRPLDE ,DRPLDT , STRAN, DSTRAN, TIME, DTIME, TEMP,
     *                DTEMP, PREDEF, DPRED, CMNAME, NDI, NSHR, NTENS,
     *                NSTATV, PROPS, NPROPS, COORDS, DROT, PNEWDT, 
     *                CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER, KSPT, 
     *                KSTEP, KINC)

C    
      INCLUDE 'ABA_PARAM.INC'
C        
	CHARACTER*80 CMNAME

	DIMENSION STRESS(NTENS), STATEV(NSTATV), DDSDDE(NTENS,NTENS),
     *     	  DDSDDT(NTENS), DRPLDE(NTENS), STRAN(NTENS), 
     *          DSTRAN(NTENS), TIME(2), PREDEF(1), DPRED(1),
     *          PROPS(NPROPS), COORDS(3), DROT(3,3), DFGRD0(3,3), 
     *          DFGRD1(3,3)
   

      integer K1,K2
      Real e0, gama, landa, PHInu, kp, ap, PHIcs, kf, kPT, aPT, Ga, Ka 
      Real Pc, q, alfa, p, pp,  Sap, ecs_p, SinPHIp, Mp, Malfa2, F, eta
      Real MPTc, SinPHIPTc, SinPHIPTe, MPTe, C, teta, Gteta, MPT, A, D 
      Real Pfin, Pcin, Rond_Pc, Rond_f, Qprime, Hn, Qzegon, R, L, dp, dq
      real Depsp, Depsq,Ecs,SA,SinPHIf, Mf,H, Pa, Pf, G, K, te,DSTRES(4)
 !!!!     
      e0 = PROPS (1)
      gama = PROPS (2)
      landa = PROPS (3)
      PHInu = PROPS (4)
      kp = PROPS (5)
      ap = PROPS (6)
      PHIcs = PROPS (7)
      kf = PROPS (8)
      kPT = PROPS (9)
      aPT = PROPS (10)
      Ga = PROPS (11)
      Ka = PROPS (12)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     
!                for undrained Depsp=0

      !!stateOld(k,1)=PC=2000000
      !!stateOld(k,2)=Mp
      !!stateOld(k,3)=pc-pf
      !!stateOld(k,4)=alfa
      !!statev(5)=counter,initial=1
      Depsp=0
      H=1.0
      Pa=101325
      e=e0
      IF (STATEV(5)==1) Then
      ![
      STATEV(5)=STATEV(5)+1
      P=(STRESS (1)+ STRESS (2)+ STRESS (3))/3+0.0001
      Depsq = 2*(DSTRAN(1)-DSTRAN(3))/3
      q= STRESS (1)- STRESS (3)
      STATEV(4)=q/STATEV(1)
      
      IF (STATEV(4)==0 .AND. (-1*DSTRAN (1))>=0) Then 
         q=1
      Else IF (STATEV(4)==0 .AND. (-1*DSTRAN (1))<0) Then
         q=-1
      End IF
      eta=q/p 
      pp = 0.64*STATEV(1)
      Ecs_p=gama-landa* LOG (pp/1000)
      SAp = e0-ecs_p
      IF ((-1*DSTRAN (1))>=0) Then 
          SinPHIp=  SIN (PHInu*3.1415/180)- kp*Sap
          STATEV(2)=(6*sinPHIp)/(3- sinPHIp)
      Else 
          SinPHIp=  SIN (PHInu*3.1415/180)- kp*Sap-ap
          STATEV(2)=(6*sinPHIp)/(3+ sinPHIp)
      End IF
      Malfa2=(5*STATEV(2)-STATEV(4))*(STATEV(2)-STATEV(4))
      F=(eta-STATEV(4))**2- Malfa2*(1-(p/STATEV(1))**0.5)
      IF (F>0) Then
          F=0 
      End IF
      Ecs =gama-landa* LOG (p/1000)
      SA = e0-ecs
      SinPHIf=  SIN (PHIcs*3.1415/180)- kf*Sa
      Mf=(6*sinPHIf)/(3- sinPHIf)
      SinPHIPTc=  SIN (PHIcs*3.1415/180)- kPT*Sa
      MPTc=(6*sinPHIPTc)/(3- sinPHIPTc)
      SinPHIPTe=  SIN (PHIcs*3.1415/180)+aPT- kPT*Sa
      MPTe=(6*sinPHIPTe)/(3- sinPHIPTe)
      !  baraye semehvari
      IF ((-1*DSTRAN (1))>=0) Then 
          MPT=MPTc
      Else 
          MPT=MPTe
      End IF
      A=9/(9-2*MPT*eta+3*MPT)
      D=(2/3)**0.5*A*( MPT-eta)
      Pf=p/(1-(Mf-STATEV(4))**2/Malfa2)**2
      G=Ga*(2.973-e)**2*(p/Pa)**0.5/(1+e)
      K=Ka*(2.973-e)**2*(p/Pa)**0.5/(1+e)
      STATEV(3)=Pf-STATEV(1)
      Rond_Pc=H*G*(Pf-STATEV(1))/STATEV(3)
      Rond_f=-1* Malfa2*(p/STATEV(1))**0.5/2/STATEV(1) 
      Qprime=2/P*(eta-STATEV(4))
      Hn=-1*(2/3)**0.5*Rond_Pc*Rond_f/Qprime
      Qzegon=(Malfa2/2/(P*STATEV(1))**0.5-2/P*(eta**2-STATEV(4)*eta))/
      R=(2/3)**0.5*3*Qzegon/Qprime
      L=(K*R*Depsp+6**0.5*G*Depsq)/(Hn+K*R*D+2*G)

      IF (L<0) Then
          L=0
      End IF
      dp=K*(Depsp-L*D)
      dq=3*G*(Depsq-(2/3)**0.5*L)
      DSTRES(1) = dp+2*dq/3
      DSTRES(2) = dp-dq/3
      DSTRES(3) = dp-dq/3
      DSTRES(4) = 0 
      DO K1=1,NTENS     
         DO K2=1,NTENS
           DDSDDE(K1, K2)=DSTRES(K1)/DSTRAN(K2)
         END DO
      End DO
      DO K1=1, NTENS     
         DO K2=1, NTENS
           STRESS (K1)= STRESS (K1)+DSTRAN(K2)* DDSDDE(K1, K2)
         END DO
      End DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      Else
  !!    [
      STATEV(3)=STATEV(3)
      STATEV(4)=STATEV(4)
      STATEV(5)=STATEV(5)+1
      P=(STRESS (1)+ STRESS (2)+ STRESS (3))/3+0.0001
      Depsq = 2*(DSTRAN(1)-DSTRAN(3))/3
      q= STRESS (1)- STRESS (3)
      if (q==0) then
            q=1
      end if
      eta=q/p
      Malfa2=(5*STATEV(2)-STATEV(4))*(STATEV(2)-STATEV(4))
      IF ((-1*DSTRAN (1))>0)  Then
          STATEV(1)=p/(1-(eta-STATEV(4))**2/STATEV(2))**2
      Else
          STATEV(1)=STATEV(1)
      End IF
      pp =STATEV(1)*(1-(STATEV(2)-STATEV(4))/(5*STATEV(2)-STATEV(4)))**2
      Ecs_p=gama-landa* LOG (pp/1000)
      SAp = e0-ecs_p
      IF ((-1*DSTRAN (1))>=0) Then 
          SinPHIp=  SIN (PHInu*3.1415/180)- kp*Sap
          STATEV(2)=(6*sinPHIp)/(3- sinPHIp)
      Else 
          SinPHIp=  SIN (PHInu*3.1415/180)- kp*Sap-ap
          STATEV(2)=(6*sinPHIp)/(3+ sinPHIp)
      End IF
      Malfa2=(5*STATEV(2)-STATEV(4))*(STATEV(2)-STATEV(4))
      F=(eta-STATEV(4))**2- Malfa2*(1-(p/STATEV(1))**0.5)
      IF (F>0) Then
          F=0 
      End IF
      Ecs =gama-landa* LOG (p/1000)
      SA = e0-ecs
      SinPHIf=  SIN (PHIcs*3.1415/180)- kf*Sa
      Mf=(6*sinPHIf)/(3- sinPHIf)
      SinPHIPTc=  SIN (PHIcs*3.1415/180)- kPT*Sa
      MPTc=(6*sinPHIPTc)/(3- sinPHIPTc)
      SinPHIPTe=  SIN (PHIcs*3.1415/180)+aPT- kPT*Sa
      MPTe=(6*sinPHIPTe)/(3- sinPHIPTe)
            !  baraye semehvari
      IF ((-1*DSTRAN (1))>=0) Then 
          MPT=MPTc
      Else 
          MPT=MPTe
      End IF

      A=9/(9-2*MPT*eta+3*MPT)
      D=(2/3)**0.5*A*(MPT-eta)
      Pf=p/(1-(Mf-STATEV(4))**2/Malfa2)**2
      G=Ga*(2.973-e)**2*(p/Pa)**0.5/(1+e)
      K=Ka*(2.973-e)**2*(p/Pa)**0.5/(1+e)
      STATEV(3)=Pf-STATEV(1)
      Rond_Pc=H*G*(Pf-STATEV(1))/STATEV(3)
      Rond_f=-1* Malfa2*(p/STATEV(1))**0.5/2/STATEV(1) 
      Qprime=2/P*(eta-STATEV(4))
      Hn=-1*(2/3)**0.5*Rond_Pc*Rond_f/Qprime
      Qzegon=(Malfa2/2/(P*STATEV(1))**0.5-2/P*(eta**2-STATEV(4)*eta))/3

      R=(2/3)**0.5*3*Qzegon/Qprime
      L=(K*R*Depsp+6**0.5*G*Depsq)/(Hn+K*R*D+2*G)

      IF (L<0) Then
          L=0
      End IF
      dp=K*(Depsp-L*D)
      dq=3*G*(Depsq-(2/3)**0.5*L)
      DSTRES(1) = dp+2*dq/3
      DSTRES(2) = dp-dq/3
      DSTRES(3) = dp-dq/3
      DSTRES(4) = 0 
      DO K1=1,NTENS     
         DO K2=1,NTENS
           DDSDDE(K1, K2)=DSTRES(K1)/DSTRAN(K2)
         END DO
      End DO
      DO K1=1, NTENS     
         DO K2=1, NTENS
           STRESS (K1)= STRESS (K1)+DSTRAN(K2)* DDSDDE(K1, K2)
         END DO
      End DO
      End IF
      return
      end