C
C User subroutine VFRIC 
      subroutine vfric (
C Write only - 
     *     fTangential, 
C Read/Write - 
     *     statev,
C Read only - 
     *     kStep, kInc, nContact, nFacNod, nSlvNod, nMstNod,
     *     nFricDir, nDir, nStateVar, nProps, nTemp, nPred, numDefTfv, 
     *     jSlvUid, jMstUid, jConSlvid, jConMstid, timStep, timGlb,
     *     dTimCur, surfInt, surfSlv, surfMast, lContType,
     *     dSlipFric, fStickForce, fTangPrev, fNormal, frictionWork,
     *     shape, coordSlv, coordMst, dircosSl, dircosN, props,
     *     areaSlv, tempSlv, preDefSlv, tempMst, preDefMst )
C
      include 'vaba_param.inc'
C
      dimension props(nProps), statev(*), 
     1     dSlipFric(nDir,nContact),
     2     fTangential(nFricDir,nContact),
     3     fTangPrev(nDir,nContact),
     4     fStickForce(nContact), areaSlv(nSlvNod),
     5     fNormal(nContact), shape(nFacNod,nContact), 
     6     coordSlv(nDir,nSlvNod), coordMst(nDir,nMstNod), 
     7     dircosSl(nDir,nContact), dircosN(nDir,nContact),
     8     jSlvUid(nSlvNod), jMstUid(nMstNod),
     9     jConSlvid(nContact), jConMstid(nFacNod,nContact),
     1     tempSlv(nContact), tempMst(numDefTfv),
     2     preDefSlv(nContact, nPred),
     3     preDefMst(numDefTfv, nPred)
C
      character*80 surfInt, surfSlv, surfMast 
      parameter ( zero = 0.d0, one = 1.d0 )
C
      xMuk = props(1)
      xMus = props(2)
      beta = props(3)
      dTimCurInv = one / dTimCur	
      if ( nDir .eq. 2 ) then
        do kcon = 1, ncontact
          fn = fNormal(kcon) 
          fs = fStickForce(kcon)
          slipRateNorm = sqrt ( 
     *         dSlipFric(1,kcon) * dSlipFric(1,kcon) +
     *         dSlipFric(2,kcon) * dSlipFric(2,kcon) ) * dTimCurInv
          xMu = xMuk + (xMus-xMuk)*exp(-beta*slipRateNorm)
          ft = min ( xMu * fn, fs ) 
          fTangential(1,kcon) = -ft
        end do
      else if ( nDir .eq. 3 ) then
        do kcon = 1, ncontact
          fn = fNormal(kcon)
          fs = fStickForce(kcon)
          slipRateNorm = sqrt ( 
     *         dSlipFric(1,kcon) * dSlipFric(1,kcon) +
     *         dSlipFric(2,kcon) * dSlipFric(2,kcon) +
     *         dSlipFric(3,kcon) * dSlipFric(3,kcon) ) * dTimCurInv
          xMu = xMuk + (xMus-xMuk)*exp(-beta*slipRateNorm)
          ft = min ( xMu * fn, fs )
          fTangential(1,kcon) = -ft
          fTangential(2,kcon) = zero
        end do
      end if
*
      return
      end