C
C User subroutine VUMAT
      subroutine vumat (
C Read only -
     *     nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     *     stepTime, totalTime, dt, cmname, coordMp, charLength,
     *     props, density, strainInc, relSpinInc,
     *     tempOld, stretchOld, defgradOld, fieldOld,
     *     stressOld, stateOld, enerInternOld, enerInelasOld,
     *     tempNew, stretchNew, defgradNew, fieldNew,
C Write only -
     *     stressNew, stateNew, enerInternNew, enerInelasNew )
C
      include 'vaba_param.inc'
C
      dimension coordMp(nblock,*), charLength(nblock), props(nprops),
     1     density(nblock), strainInc(nblock,ndir+nshr),
     2     relSpinInc(nblock,nshr), tempOld(nblock),
     3     stretchOld(nblock,ndir+nshr), 
     4     defgradOld(nblock,ndir+nshr+nshr),
     5     fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
     6     stateOld(nblock,nstatev), enerInternOld(nblock),
     7     enerInelasOld(nblock), tempNew(nblock),
     8     stretchNew(nblock,ndir+nshr),
     9     defgradNew(nblock,ndir+nshr+nshr),
     1     fieldNew(nblock,nfieldv),
     2     stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     3     enerInternNew(nblock), enerInelasNew(nblock)
C     
      character*80 cmname
      parameter ( zero = 0.d0, one = 1.d0, two = 2.d0,
     *     third = 1.d0 / 3.d0, half = 0.5d0, op5 = 1.5d0 )
C
C For plane strain, axisymmetric, and 3D cases using
C the J2 Mises Plasticity with linear hardening.
C The state variable is stored as:
C            STATE(*,1) = equivalent plastic strain
C            STATE(*,2) = yield stress
C
C User needs to input
C     props(1)      Young's modulus
C     props(2)      Poisson's ratio
C     props(3)      Yield stress for temperatures <= T1
C     props(4)      Temperature T1
C     props(5)      Yield stress for temperatures >= T2
C     props(6)      Temperature T2
C
      e      = props(1)
      xnu    = props(2)
      yield1 = props(3)
      temp1  = props(4)
      yield2 = props(5)
      temp2  = props(6)
      tempSoft = (yield2 - yield1)/(temp2 - temp1)  
      twomu  = e / ( one + xnu )
      alamda = xnu * twomu / ( one - two * xnu )
      thremu = op5 * twomu
*
      if ( stepTime .eq. zero ) then    
        do k = 1, nblock
          trace = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
          stressNew(k,1) = stressOld(k,1) 
     *         + twomu * strainInc(k,1) + alamda * trace
          stressNew(k,2) = stressOld(k,2) 
     *         + twomu * strainInc(k,2) + alamda * trace
          stressNew(k,3) = stressOld(k,3) 
     *         + twomu * strainInc(k,3) + alamda * trace
          stressNew(k,4)=stressOld(k,4) + twomu * strainInc(k,4)
          if ( nshr .gt. 1 ) then
            stressNew(k,5)=stressOld(k,5) + twomu * strainInc(k,5)
            stressNew(k,6)=stressOld(k,6) + twomu * strainInc(k,6)
          end if
        end do
      else
        do k = 1, nblock
          yieldNew = zero
          if ( tempNew(k) .le. temp1 ) then
             yieldNew = yield1
          else if ( tempNew(k) .ge. temp2 ) then
             yieldNew = yield2
          else
             yieldNew = yield1 + tempSoft * (tempNew(k)-temp1)
          end if          
          trace = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
          s11 = stressOld(k,1) + twomu * strainInc(k,1) + alamda * trace
          s22 = stressOld(k,2) + twomu * strainInc(k,2) + alamda * trace
          s33 = stressOld(k,3) + twomu * strainInc(k,3) + alamda * trace 
          s12 = stressOld(k,4) + twomu * strainInc(k,4)
          if ( nshr .gt. 1 ) then
            s13 = stressOld(k,5) + twomu * strainInc(k,5)
            s23 = stressOld(k,6) + twomu * strainInc(k,6)
          end if
*     
          smean = third * ( s11 + s22 + s33 ) 
          s11 = s11 - smean
          s22 = s22 - smean
          s33 = s33 - smean
          if ( nshr .eq. 1 ) then
            vmises = sqrt( op5*(s11*s11+s22*s22+s33*s33+two*s12*s12) )
          else
            vmises = sqrt( op5 * ( s11 * s11 + s22 * s22 + s33 * s33 +
     *           two * s12 * s12 + two * s13 * s13 + two * s23 * s23 ) )
          end if
*
          sigdif = vmises - yieldNew
          facyld = zero
          if ( sigdif .gt. zero ) facyld = one
          deqps = facyld * sigdif / thremu 
*
* Update the state variable
          stateNew(k,1) = stateOld(k,1) + deqps
          stateNew(k,2) = yieldNew
*
* Update the stress
          factor = yieldNew / ( yieldNew + thremu * deqps )
          stressNew(k,1) = s11 * factor + smean
          stressNew(k,2) = s22 * factor + smean
          stressNew(k,3) = s33 * factor + smean
          stressNew(k,4) = s12 * factor
          if ( nshr .gt. 1 ) then
            stressNew(k,5) = s13 * factor
            stressNew(k,6) = s23 * factor
          end if
*     
* Update the specific internal energy -
          if ( nshr .eq. 1 ) then
            stressPower = half * (
     *        ( stressOld(k,1) + stressNew(k,1) ) * strainInc(k,1) +
     *        ( stressOld(k,2) + stressNew(k,2) ) * strainInc(k,2) +
     *        ( stressOld(k,3) + stressNew(k,3) ) * strainInc(k,3) ) +
     *        ( stressOld(k,4) + stressNew(k,4) ) * strainInc(k,4)     
          else
            stressPower = half * (
     *        ( stressOld(k,1) + stressNew(k,1) ) * strainInc(k,1) +
     *        ( stressOld(k,2) + stressNew(k,2) ) * strainInc(k,2) +
     *        ( stressOld(k,3) + stressNew(k,3) ) * strainInc(k,3) ) +
     *        ( stressOld(k,4) + stressNew(k,4) ) * strainInc(k,4) +    
     *        ( stressOld(k,5) + stressNew(k,5) ) * strainInc(k,5) +    
     *        ( stressOld(k,6) + stressNew(k,6) ) * strainInc(k,6)     
          end if
          enerInternNew(k) = enerInternOld(k) + stressPower / density(k)
*     
* Update the dissipated inelastic specific energy -
          yieldOld = stateOld(k,2)
          plasticWorkInc = half * ( yieldOld + yieldNew ) * deqps
          enerInelasNew(k) = enerInelasOld(k)
     *         + plasticWorkInc / density(k)
        end do
      end if
*     
      return
      end