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

	  INCLUDE 'ABA_PARAM.INC'

      CHARACTER*80 CMNAME
      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),DFGRD0(3,3),DFGRD1(3,3)

      INTEGER NProny, counter
      REAL*8 rhoi(7)
      REAL*8 C11inf, C12inf, C22inf, C23inf, C44inf, C55inf 
      REAL*8 C11i(7), C12i(7), C22i(7), C23i(7), C44i(7), C55i(7)
      REAL*8 qold(NTENS, 7), q(NTENS, 7)
      REAL*8 f111, f122, f123, f121, f222, f233, f232, f223
      REAL*8 f444, f555, f556
      REAL*8 depsilon(NTENS), epsilonE(NTENS), sigmaE(NTENS)    
          
C----------PROPS DESIGNATION----------C
C
C     PROPS  1 - 7 : rhoi   %relaxation times
C     PROPS  8      : C11inf %
C     PROPS  9 - 15 : C11i   %       
C     PROPS 16      : C12inf %
C     PROPS 17 - 23 : C12i   %
C     PROPS 24      : C22inf %
C	PROPS 25 - 31 : C22i   %
C     PROPS 32      : C23inf %
C     PROPS 33 - 39 : C23i   %
C     PROPS 40      : C44inf %
C     PROPS 41 - 47 : C44i   %
C     PROPS 48      : C55inf %
C     PROPS 49 - 55 : C55i   %

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

C	Number of Prony terms
      NProny=7 
            
C     Obtain Prony coefficients 
      do i=1,NProny
        rhoi(i)=PROPS(i)
        C11i(i)=PROPS(8+i)
        C12i(i)=PROPS(16+i)
        C22i(i)=PROPS(24+i)
        C23i(i)=PROPS(32+i)
        C44i(i)=PROPS(40+i)
        C55i(i)=PROPS(48+i)
      end do
	  
C	Obtain long term modulus
	C11inf=PROPS(8)
	C12inf=PROPS(16)
	C22inf=PROPS(24)
	C23inf=PROPS(32)
	C44inf=PROPS(40)
	C55inf=PROPS(48)
	  
C----------solution-dependent variables----------C
	  
C     old heredity integrals
      counter = 1 
	do j = 1,NTENS
      	do i = 1,NProny
	    qold(j,i) = STATEV(counter) 
	    counter = counter+1              
	  end do
      end do
	        	       	  
C----------update stress--------------C
      
C     switch to Voight notation      
      depsilon(1)=DSTRAN(1)
      depsilon(2)=DSTRAN(2)
      depsilon(3)=DSTRAN(3)
      depsilon(4)=DSTRAN(6)
      depsilon(5)=DSTRAN(5)
      depsilon(6)=DSTRAN(4)
      epsilonE(1)=STRAN(1) + DSTRAN(1)
      epsilonE(2)=STRAN(2) + DSTRAN(2)
      epsilonE(3)=STRAN(3) + DSTRAN(3)
      epsilonE(4)=STRAN(6) + DSTRAN(6)
      epsilonE(5)=STRAN(5) + DSTRAN(5)
      epsilonE(6)=STRAN(4) + DSTRAN(4)
            
C     compute q
      do j = 1,NTENS
        do i = 1, NProny	
          q(j,i)=EXP(-DTIME/rhoi(i))*qold(j,i)+ 
     & depsilon(j)*(1-EXP(-DTIME/rhoi(i)))/(DTIME/rhoi(i))
        end do
      end do
      
C     compute f
	f111=0
      f122=0
      f123=0
      f121=0
      f222=0
      f233=0
      f232=0
      f223=0
      f444=0
      f555=0
      f556=0
      do i = 1, NProny
        f111=f111+C11i(i)*q(1,i)
        f122=f122+C12i(i)*q(2,i)
        f123=f123+C12i(i)*q(3,i)
        f121=f121+C12i(i)*q(1,i)
        f222=f222+C22i(i)*q(2,i)
        f233=f233+C23i(i)*q(3,i)
        f232=f232+C23i(i)*q(2,i)
        f223=f223+C22i(i)*q(3,i)
        f444=f444+C44i(i)*q(4,i)
        f555=f555+C55i(i)*q(5,i)
        f556=f556+C55i(i)*q(6,i)
      end do
      
      sigmaE(1) = C11inf*epsilonE(1)+f111+C12inf*epsilonE(2)+f122+ 
     & C12inf*epsilonE(3)+f123   
      sigmaE(2) = C12inf*epsilonE(1)+f121+C22inf*epsilonE(2)+f222+ 
     & C23inf*epsilonE(3)+f233
      sigmaE(3) = C12inf*epsilonE(1)+f121+C23inf*epsilonE(2)+f232+ 
     & C22inf*epsilonE(3)+f223
      sigmaE(4) = C44inf*epsilonE(4) + f444
      sigmaE(5) = C55inf*epsilonE(5) + f555
      sigmaE(6) = C55inf*epsilonE(6) + f556      
            
      STRESS(1) = sigmaE(1)
      STRESS(2) = sigmaE(2)
      STRESS(3) = sigmaE(3)
      STRESS(4) = sigmaE(6)
      STRESS(5) = sigmaE(5)
      STRESS(6) = sigmaE(4)

C-----update solution-dependent variables-----C

C     hereditary integrals 
      counter = 1
      do j = 1,NTENS
        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, NTENS
        do j = 1, NTENS
          DDSDDE(i,j)=0.0  
        end do
      end do

      do i = 1, NProny
        DDSDDE(1,1)=DDSDDE(1,1)+C11inf+
     & C11i(i)*(1-EXP(-DTIME/rhoi(i)))/(DTIME/rhoi(i))
        DDSDDE(1,2)=DDSDDE(1,2)+C12inf+
     & C12i(i)*(1-EXP(-DTIME/rhoi(i)))/(DTIME/rhoi(i))
        DDSDDE(2,2)=DDSDDE(2,2)+C22inf+
     & C22i(i)*(1-EXP(-DTIME/rhoi(i)))/(DTIME/rhoi(i))
        DDSDDE(2,3)=DDSDDE(2,3)+C23inf+
     & C23i(i)*(1-EXP(-DTIME/rhoi(i)))/(DTIME/rhoi(i))
        DDSDDE(5,5)=DDSDDE(5,5)+C55inf+
     & C55i(i)*(1-EXP(-DTIME/rhoi(i)))/(DTIME/rhoi(i))
        DDSDDE(6,6)=DDSDDE(6,6)+C44inf+
     & C44i(i)*(1-EXP(-DTIME/rhoi(i)))/(DTIME/rhoi(i))
      end do

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

C     update elastic strain energy in SSE
      SSE = 0.0
    
C     update creep dissipation in SCD
      SCD = 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