!
! User subroutine to define the frictional coefficient for contact surfaces.
!
      subroutine fric_coef(fCoef, fCoefDeriv, 
     *   nBlock, nProps, nTemp, nFields, 
     *   jFlags, rData, surfInt, surfSec, surfMain, 
     *   props, slipData, pressure, tempAvg, fieldAvg )

      include 'aba_param.inc'

      real(nprecd*4), intent(out) :: fCoef(nBlock)         ! friction coefficient at the contact point
      real(nprecd*4), intent(out) :: fCoefDeriv(nBlock,20)  ! derivatives of the friction coefficient

      integer, intent(in)         :: nBlock, nProps, nTemp, nFields

      real(nprecd*4), parameter   :: ZERO = 0.0, ONE = 1.0

      parameter( jKStep   = 1,
     *           jKInc    = 2,
     *           nFlags   = 2,
     *           jTimStep = 1,
     *           jTimGlb  = 2,
     *           jDTimCur = 3,
     *           nData    = 3,
     *           idfdSlipRate = 1,
     *           idfdPressure = 2,
     *           idfdtemperature = 3)

      dimension  props(nProps),slipData(nBlock,6),pressure(nBlock),
     *   tempAvg(nBlock), fieldAvg(nBlock,nFields),
     *   jFlags(nFlags), rData(nData)

      character*80 surfInt, surfSec, surfMain 
      
      call aset(fCoef,ZERO,nBlock)
      call aset(fCoefDeriv,ZERO, nBlock*20)

      Ck = props(1)
      Dk = props(2)
      Cs = props(3)
      Ds = props(4)
      dc = props(5)     
      do k1=1,nBlock
        uk    = Ck*tempAvg(k1)+Dk*pressure(k1)        
        dukdt = Ck 
        dukdp = Dk
        us    = Cs*tempAvg(k1)+Ds*pressure(k1)        
        dusdt = Cs 
        dusdp = Ds
        decay = exp(-dc*slipData(k1,1))
        u     = uk+(us-uk)*decay
        duduk = ONE-decay
        dudus = decay   
        fCoef(k1) = u        
        fCoefDeriv(k1,idfdSlipRate) = -dc*(us-uk)*decay
        fCoefDeriv(k1,idfdPressure) = duduk*dukdp+dudus*dusdp
        fCoefDeriv(k1,idfdTemperature) = duduk*dukdt+dudus*dusdt
      end do

      return
      end