*USER SUBROUTINES
C***********************************************************************
C23456789012345678901234567890123456789012345678901234567890123456789012
      SUBROUTINE UGENS(DDNDDE,FORCE,STATEV,SSE,SPD,PNEWDT,STRAN,
     1 DSTRAN,TSS,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CENAME,NDI,
     2 NSHR,NSECV,NSTATV,PROPS,JPROPS,NPROPS,NJPROP,COORDS,CELENT,
     3 THICK,DFGRD,CURV,BASIS,NOEL,NPT,KSTEP,KINC,NIT,LINPER)

      INCLUDE 'ABA_PARAM.INC'

      CHARACTER*80 CMNAME
      DIMENSION DDNDDE(NSECV,NSECV),FORCE(NSECV),STATEV(NSTATV),
     1 STRAN(NSECV),DSTRAN(NSECV),TSS(2),TIME(2),PREDEF(*),
     2 DPRED(*),PROPS(*),JPROPS(*),COORDS(3),DFGRD(3,3),
     3 CURV(2,2),BASIS(3,3)

      INTEGER nProny, counter
      REAL*8 rhoi(7)
      REAL*8 A11inf, A12inf, A33inf, D11inf, D12inf, D33inf 
      REAL*8 A11i(7), A12i(7), A33i(7), D11i(7), D12i(7), D33i(7) 
      REAL*8 qold(NSECV, 7), q(NSECV, 7)
      REAL*8 a111, a122, a121, a112, a333
      REAL*8 d114, d125, d124, d115, d336
      REAL*8 strainE(NSECV)
      REAL*8 c1, c2, Tr, aT
          
C----------PROPS DESIGNATION----------C
C
C     PROPS  1 - 7  : rhoi   %relaxation times
C     PROPS  8      : A11inf %
C     PROPS  9 - 15 : A11i   %       
C     PROPS 16      : A12inf %
C     PROPS 17 - 23 : A12i   %
C     PROPS 24      : A33inf %
C	PROPS 25 - 31 : A33i   %
C     PROPS 32      : D11inf %
C     PROPS 33 - 39 : D11i   %
C     PROPS 40      : D12inf %
C     PROPS 41 - 47 : D12i   %
C     PROPS 48      : D33inf %
C     PROPS 49 - 55 : D33i   %

C----------STATEV DESIGNATION----------C
C
C     STATEV  1 - 42 : qold   %old hereditary integrals
	      
C-----material parameters---------C

C	Number of Prony terms
      nProny=7 
            
C     Obtain Prony coefficients 
      do i=1,nProny
        rhoi(i)=PROPS(i)
        A11i(i)=PROPS(8+i)
        A12i(i)=PROPS(16+i)
        A33i(i)=PROPS(24+i)
        D11i(i)=PROPS(32+i)
        D12i(i)=PROPS(40+i)
        D33i(i)=PROPS(48+i)
      end do
	  
C	Obtain long term modulus
	A11inf=PROPS(8)
	A12inf=PROPS(16)
	A33inf=PROPS(24)
	D11inf=PROPS(32)
	D12inf=PROPS(40)
	D33inf=PROPS(48)
	  
C	WLF parameters
      Tr=313.0
      c1=28.3816
	c2=93.291
      
C----------solution-dependent variables----------C
	  
C     old heredity integrals
      counter = 1 
	do j = 1,NSECV
      	do i = 1,nProny
	    qold(j,i) = STATEV(counter) 
	    counter = counter+1              
	  end do
      end do
	        	       	  
C----------update stress--------------C
      
C     generalized section strains at the end of increment     
      strainE(1)=STRAN(1) + DSTRAN(1)
      strainE(2)=STRAN(2) + DSTRAN(2)
      strainE(3)=STRAN(3) + DSTRAN(3)
      strainE(4)=STRAN(4) + DSTRAN(4)
      strainE(5)=STRAN(5) + DSTRAN(5)
      strainE(6)=STRAN(6) + DSTRAN(6)
      
C	shift factor    
      aT=10**(-c1*(TEMP-Tr)/(c2+TEMP-Tr))
                  
C     reduced time increment
	dtr=DTIME/aT
      
C     compute q
      do j = 1,NSECV
        do i = 1, nProny	
          q(j,i)=exp(-dtr/rhoi(i))*qold(j,i)+ 
     & DSTRAN(j)*(1-exp(-dtr/rhoi(i)))/(dtr/rhoi(i))
        end do
      end do
      
C     compute a and d
 	a111=0
      a122=0
      a121=0
      a112=0
      a333=0
      d114=0
      d125=0
      d124=0
      d115=0
      d336=0
      do i = 1, nProny
        a111=a111+A11i(i)*q(1,i)
        a122=a122+A12i(i)*q(2,i)
        a121=a121+A12i(i)*q(1,i)
        a112=a112+A11i(i)*q(2,i)
        a333=a333+A33i(i)*q(3,i)
        d114=d114+D11i(i)*q(4,i)
        d125=d125+D12i(i)*q(5,i)
        d124=d124+D12i(i)*q(4,i)
        d115=d115+D11i(i)*q(5,i)
        d336=d336+D33i(i)*q(6,i)
      end do
            
      FORCE(1) = A11inf*strainE(1)+a111+A12inf*strainE(2)+a122 
      FORCE(2) = A12inf*strainE(1)+a121+A11inf*strainE(2)+a112 
      FORCE(3) = A33inf*strainE(3)+a333
      FORCE(4) = D11inf*strainE(4)+d114+D12inf*strainE(5)+d125
C      FORCE(4) = D11inf*strainE(4)+d114
      FORCE(5) = D12inf*strainE(4)+d124+D11inf*strainE(5)+d115
C      FORCE(5) = D11inf*strainE(5)+d115
      FORCE(6) = D33inf*strainE(6)+d336
                  
C-----update solution-dependent variables-----C

C     hereditary integrals 
      counter = 1
      do j = 1,NSECV
        do i = 1, nProny	
          STATEV(counter)=q(j,i)
          counter=counter+1
        end do
      end do

C-----update Jacobian (tangent stiffess)-----C

      do i = 1, NSECV
        do j = 1, NSECV
          DDNDDE(i,j)=0.0  
        end do
      end do

      do i = 1, nProny
        DDNDDE(1,1)=DDNDDE(1,1)+A11inf+
     & A11i(i)*(1-exp(-dtr/rhoi(i)))/(dtr/rhoi(i))
        DDNDDE(1,2)=DDNDDE(1,2)+A12inf+
     & A12i(i)*(1-exp(-dtr/rhoi(i)))/(dtr/rhoi(i))
        DDNDDE(3,3)=DDNDDE(3,3)+A33inf+
     & A33i(i)*(1-exp(-dtr/rhoi(i)))/(dtr/rhoi(i))
        DDNDDE(4,4)=DDNDDE(4,4)+D11inf+
     & D11i(i)*(1-exp(-dtr/rhoi(i)))/(dtr/rhoi(i))
        DDNDDE(4,5)=DDNDDE(4,5)+D12inf+
     & D12i(i)*(1-exp(-dtr/rhoi(i)))/(dtr/rhoi(i))
        DDNDDE(6,6)=DDNDDE(6,6)+D33inf+
     & D33i(i)*(1-exp(-dtr/rhoi(i)))/(dtr/rhoi(i))
      end do

      DDNDDE(2,1)=DDNDDE(1,2)
      DDNDDE(2,2)=DDNDDE(1,1)
      DDNDDE(5,4)=DDNDDE(4,5)
      DDNDDE(5,5)=DDNDDE(4,4)
            
C----------update energy----------C

C     update elastic strain energy in SSE
      SSE = 0.0
    
C     update plastic dissipation in SPD
	SPD = 0.0
	
      RETURN
      END
C***********************************************************************
C***********************************************************************
C23456789012345678901234567890123456789012345678901234567890123456789012
      SUBROUTINE SDVINI(STATEV,COORDS,NSTATV,NCRDS,NOEL,NPT,
     1 LAYER,KSPT)

      INCLUDE 'ABA_PARAM.INC'

      DIMENSION STATEV(NSTATV),COORDS(NCRDS) 

C     Initialize state variables
      do i = 1,NSTATV
        STATEV(i) = 0.0
      end do

      RETURN
      END