c
c User subroutine VUSDFLD for user-defined fields
c
      subroutine vusdfld(
c Read only -
     *   nblock, nstatev, nfieldv, nprops, ndir, nshr, 
     *   jElemUid, kIntPt, kLayer, kSecPt, 
     *   stepTime, totalTime, dt, cmname, 
     *   coordMp, direct, T, charLength, props, 
     *   stateOld, 
c Write only -
     *   stateNew, field )
c
      include 'vaba_param.inc'
c
      dimension props(nprops),
     *          jElemUid(nblock), coordMp(nblock, *), 
     *          direct(nblock, 3, 3), T(nblock,3,3), 
     *          charLength(nblock),
     *          stateOld(nblock, nstatev), 
     *          stateNew(nblock, nstatev),
     *          field(nblock, nfieldv)
      character*80 cmname
c Properties array
c     props(1) -> eqpsCrit
c     props(2) -> failDisp
c     props(3) -> Dmax
c
      character*3 cData(maxblk)
      dimension jData(maxblk)
      dimension eqps(maxblk)
c
      parameter ( zero = 0.d0 )
c
c Multi-process information (if needed)
      call vgetnumcpus( numProcesses )
      call vgetrank( kProcessNum )
c
c Read properties
      eqpsCrit = props(1)
      failDisp = props(2)
      Dmax = props(3)
c Get stresses and strains from previous increment
      jStatus = 1
      call vgetvrm( 'PEEQ', eqps, jData, cData, jStatus )
c
      do k = 1, nblock
         eqpsOld = stateOld(k,1)
         damageOld = stateOld(k,2)
         damageNew = damageOld
         if ( eqps(k) .gt. eqpsCrit ) then
            deqps = eqps(k) - eqpsOld
            if ( damageOld .eq. zero ) ! First time criterion is met
     *           deqps = eqps(k) - eqpsCrit
            dUeq = charLength(k) * deqps 
            damageNew = damageOld + dUeq / failDisp
         end if
         if ( damageNew .ge. Dmax ) then
c     Element deletion
            stateNew(k,3) = zero
            damageNew = Dmax
         end if
         stateNew(k,1) = eqps(k)
         stateNew(k,2) = damageNew
         field(k,1) = damageNew
      end do
c
      return
      end