C---- SUB. CONTACT_FORCE

      SUBROUTINE CONTACT_FORCE

     &          (TIME,UPAR,NPAR,PEN,RVEL,JFLAG,IFLAG,RESULT)

C---- TO EXPORT * SUBROUTINE

      !DEC$ ATTRIBUTES DLLEXPORT,C::CONTACT_FORCE

C---- INCLUDE SYSTEM CALL

      INCLUDE 'SYSCAL.F'

C---- DEFINE VARIABLES

C     Parameter Information

C     TIME   : Simulation time of RD/Solver. (Input)

C     UPAR   : Parameters defined by user. (Input)

C     NPAR   : Number of user parameters. (Input)

C     PEN    : Penetration at the contact point. (Input)

C     RVEL   : Relative velocity between contact pair w.r.t. contact reference frame. (Input, Size : 3)

C     JFLAG  : When RD/Solver evaluates a Jacobian, the flag is true. (Input)

C     IFLAG  : When RD/Solver initializes arrays, the flag is true. (Input)

C     RESULT : Returned contact force vector. (Output, Size : 6)

      DOUBLE PRECISION TIME, UPAR(*), PEN, RVEL(*)

      INTEGER NPAR

      LOGICAL JFLAG, IFLAG

      DOUBLE PRECISION RESULT[REFERENCE](6)

C---- USER STATEMENT

      DOUBLE PRECISION DPEN, K, C, M1, M2, M3, STV, DTV, SFC, DFC

      DOUBLE PRECISION SFO, DFO, TVEL, FC(3), ZERO, MU

      DOUBLE PRECISION FRICTIONFORCE

      INTEGER I

      LOGICAL ERRFLG

      DPEN = RVEL(3)

      K = UPAR(1)

      C = UPAR(2)

      M1 = UPAR(3)

      M2 = UPAR(4)

      M3 = UPAR(5)

      STV = UPAR(6)

      DTV = UPAR(7)

      SFC = UPAR(8)

      DFC = UPAR(9)

      ZERO = 0.0D0

C---- COMPUTE CONTACT NORMAL FORCE

      IF ( DABS(PEN) .LT. 1.0D-12 ) THEN

         SFO = 0.0D0

      ELSE

         SFO = K * (DABS(PEN)**M1)

      ENDIF

      IF ( DABS(DPEN) .LT. 1.0D-12 ) THEN

         DFO = 0.0D0

      ELSE

         DFO = C * (DABS(PEN)**M3) * (DABS(DPEN)**M2)

      ENDIF

      IF ( DPEN .GT. ZERO ) DFO = - DFO

      FC(3) = SFO + DFO

C---- COMPUTE FRICTION FORCE

      TVEL = SQRT(RVEL(1)*RVEL(1) + RVEL(2)*RVEL(2))

      FC(1) = 0.0D0

      FC(2) = 0.0D0

      IF ( TVEL .GE. 1.0D-12 ) THEN

         IF ( TVEL .LT. STV ) THEN

            CALL RD_HAVSIN(DABS(TVEL),-STV,-SFC,STV,SFC,0,MU,ERRFLG)

         ELSE

            CALL RD_HAVSIN(DABS(TVEL),STV,SFC,DTV,DFC,0,MU,ERRFLG)

         ENDIF

         FRICTIONFORCE = MU * DABS(FC(3))

         IF ( TVEL == 0.0d0 ) THEN

            FC(1) = 0.0D0

            FC(2) = 0.0D0

         ELSE

            FC(1) = -FRICTIONFORCE*RVEL(1)/TVEL

            FC(2) = -FRICTIONFORCE*RVEL(2)/TVEL

         ENDIF

      ENDIF

      DO I = 1, 3

        RESULT(I) = FC(I)

      ENDDO

      RETURN

      END