SUBROUTINE SDVINI(STATEV,COORDS,NSTATV,NCRDS,NOEL,NPT,LAYER,KSPT)
	  INCLUDE 'ABA_PARAM.INC'
	  DIMENSION STATEV(NSTATV),COORDS(NCRDS)
	  
	  integer i
      
	  do i=1,NSTATV
         STATEV(i) = 0
      end do
	  
	  if (NSTATV == 4) then
	    STATEV(4) = 1 ! 1 - material point is active, 0 - is to be deleted. flag of element deletion, if it is defined
	  endif
      
	  RETURN
      END SUBROUTINE SDVINI
      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,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)

      INCLUDE 'ABA_PARAM.INC'

      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)
     
c*********************************************************************
c     Progressive Failure model for sand
c
c     Version: 8 || applicable for 3D stress state only
c
c     ----- Описание -----
c
c     Массив PROPS:
c
c     (1)     K(1)
c     (2)     G(1)
c     (3)     C_0
c     (4)     Phi_0
c
c     Массив STATEV:
c
c     (1)     Failure Flag
c     (2)     Failure Index
c     (3)                  приращения деформаций немеханической природы
c     (4)     Failure Flag element deletion (1 is actived 0 is deleted)
c
c*********************************************************************

      
       
          
      
      
      integer             i, suggest_bisection
      
      double precision, allocatable :: Total_Load(:)
	  
	  double precision, allocatable :: K(:), G(:)
	  !double precision    K(1), G(1)
	  
	  double precision, allocatable :: Max_Failure_Inc(:), Deg_Factor(:)
	  !double precision    Max_Failure_Inc, Deg_Factor
	  
	  double precision, allocatable :: eps(:), eps_kk(:)
	  !double precision    eps (NTENS), eps_kk
	  
	  double precision, allocatable :: Csh(:), Phi(:)
	  !double precision    C_0, Phi_0
	  
	  double precision, allocatable :: SP(:)
	  !double precision    s1, s3
	  
	  double precision, allocatable :: EVALS(:)
      !double precision, dimension(3) :: EVALS
c**************************************
c  Настройки, задаваемые пользователем:
c**************************************
c Максимальный инкремент шага по времени, когда происходит разрушение
      allocate(Max_Failure_Inc(1))
	  Max_Failure_Inc(1) = 0.001
c Множитель деградации модуля сдвига в случае разрушения
      allocate(Deg_Factor(1))
	  Deg_Factor(1) = 0.01
c Величина полного давления
      allocate(Total_Load(1))
	  Total_Load(1) = 200000
      
c**************************************
c     Установка флага рекомендации бисекции в положение 0 (не нужна)
c**************************************
      
      suggest_bisection = 0
      
c**************************************
c     Обнуление приращения деформаций немеханической природы
c**************************************
      
      STATEV(3) = 0
      
c**************************************
c     Задание материальных констант
c**************************************
      allocate(K(1))
	  allocate(G(1))
      K(1) = PROPS(1)
      G(1) = PROPS(2)
      
c**************************************
c     Задание констант критерия разрушения
c**************************************
      allocate(Csh(1))
	  allocate(Phi(1))     
      Csh(1) = PROPS(3)
      Phi(1) = PROPS(4)

c**************************************
c     Вычисление актуальных деформаций
c**************************************
      allocate(eps(NTENS))
      eps = STRAN + DSTRAN
      
c**************************************
c     Корректировка G(1), если уже есть разрушение
c**************************************

      if (STATEV(1) == 1) then
          G(1) = G(1)*Deg_Factor(1)
      endif
      
c**************************************
c     Первоначальное вычисление напряжений
c**************************************
      allocate(eps_kk(1))
      eps_kk(1) = eps(1) + eps(2) + eps(3)
      
      STRESS(1) = (K(1) - G(1)*2.0/3.0)*eps_kk(1) + 2.0*G(1)*eps(1)
      STRESS(2) = (K(1) - G(1)*2.0/3.0)*eps_kk(1) + 2.0*G(1)*eps(2)
      STRESS(3) = (K(1) - G(1)*2.0/3.0)*eps_kk(1) + 2.0*G(1)*eps(3)
      
      STRESS(4) = G(1)*eps(4)
      STRESS(5) = G(1)*eps(5)
      STRESS(6) = G(1)*eps(6)
      
c**************************************
c     Вычисление главных напряжений
c**************************************
	  	
	  allocate(EVALS(3))	  
	  CALL SPRINC(STRESS,EVALS,1,NDI,NSHR)
	  
	  allocate(SP(3))
      SP(1) = max(EVALS(1), EVALS(2), EVALS(3))
	  SP(2) = 0
      SP(3) = min(EVALS(1), EVALS(2), EVALS(3))
     
      deallocate(EVALS)	 
c**************************************
c     Индекс разрушения
c**************************************

      if (STATEV(1) == 0) then
!          STATEV(2) = ((s1-s3)/2*cosd(Phi_0) +
!     &     ((s1+s3)/2 + (s1-s3)/2*sind(Phi_0))*tand(Phi_0))/C_0
          STATEV(2) = 1/(2.0*Csh(1))*( (SP(1)-SP(3))/cosd(Phi(1)) +
     &     (SP(1)+SP(3))*tand(Phi(1)) )
      endif
      deallocate(Csh)
	  deallocate(Phi)
      deallocate(SP)
c**************************************
c     Проверка критерия разрушения
c**************************************

      if ((STATEV(2) >= 1) .and. (STATEV(1) == 0)) then
      
!         Флаг рекомендации бисекции установлен в положение 1 (желательна)
          suggest_bisection = 1

!         Флаг разрушения
          STATEV(1) = 1
!         Индекс разрушения
          STATEV(2) = 1
!         Флаг удаления элемента (удаление элемента = 0 )
		  STATEV(4) = 0
!         Деградация сдвиговых свойств
          G(1) = G(1)*Deg_Factor(1)
          
c         Пересчет напряжений
          STRESS(1) = (K(1) - G(1)*2.0/3.0)*eps_kk(1) + 2.0*G(1)*eps(1)
          STRESS(2) = (K(1) - G(1)*2.0/3.0)*eps_kk(1) + 2.0*G(1)*eps(2)
          STRESS(3) = (K(1) - G(1)*2.0/3.0)*eps_kk(1) + 2.0*G(1)*eps(3)
          STRESS(4) = G(1)*eps(4)
          STRESS(5) = G(1)*eps(5)
          STRESS(6) = G(1)*eps(6)
          
c         Учет деформаций немеханической природы, которые произошли с момента начала нагружения
         STATEV(3) = TIME(2)*Total_Load(1)/K(1)
      
      endif
	  deallocate(Deg_Factor)
      deallocate(eps_kk)
	  deallocate(eps)
c**************************************
c     Вычисление приращения деформаций немеханической природы
c**************************************
      !STATEV(3) = DTIME*Total_Load(1)/K(1)
      if (STATEV(1) == 1) then
      
         !STATEV(3) = STATEV(3) + DTIME*Total_Load(1)/K(1)
		STATEV(3) = DTIME*Total_Load(1)/K(1)
      
      endif
      deallocate(Total_Load)
c**************************************
c     Вычисление якобиана
c**************************************

      DDSDDE(1,1) = (K(1) - G(1)*2.0/3.0) + 2.0*G(1)
      DDSDDE(2,2) = (K(1) - G(1)*2.0/3.0) + 2.0*G(1)
      DDSDDE(3,3) = (K(1) - G(1)*2.0/3.0) + 2.0*G(1)
      DDSDDE(4,4) = G(1)
      DDSDDE(5,5) = G(1)
      DDSDDE(6,6) = G(1)
      DDSDDE(1,2) = (K(1) - G(1)*2.0/3.0)
      DDSDDE(1,3) = (K(1) - G(1)*2.0/3.0)
      DDSDDE(2,3) = (K(1) - G(1)*2.0/3.0)
      DDSDDE(2,1) = (K(1) - G(1)*2.0/3.0)
      DDSDDE(3,1) = (K(1) - G(1)*2.0/3.0)
      DDSDDE(3,2) = (K(1) - G(1)*2.0/3.0)
      deallocate(K)
	  deallocate(G)

c**********************************************
c     Обновление удельной энергии
c     упругой деформации
c**********************************************
        
      do i=1,NTENS
          SSE = SSE + STRESS(i)*DSTRAN(i)
      enddo

c**************************************
c     Управление шагом решения (рекомендация Абакусу)
c**************************************

      if ((suggest_bisection == 1) .and. (DTIME > Max_Failure_Inc(1))) then
          PNEWDT=0.5
      endif
      deallocate(Max_Failure_Inc) 
      
	  
      RETURN

      END SUBROUTINE UMAT

c**************************************
c
c**************************************
      
!      SUBROUTINE UEXPAN(EXPAN,DEXPANDT,TEMP,TIME,DTIME,PREDEF,
!     & DPRED,STATEV,CMNAME,NSTATV,NOEL)

!      INCLUDE 'ABA_PARAM.INC'

!      CHARACTER*80 CMNAME

!      DIMENSION EXPAN(*),DEXPANDT(*),TEMP(2),TIME(2),PREDEF(*),
!     & DPRED(*),STATEV(NSTATV)
      
!      EXPAN(1) = 0!STATEV(3)
	  
!      RETURN
!      END SUBROUTINE UEXPAN