C---- SUB. MATRIX_FORCE

      SUBROUTINE MATRIX_FORCE

     &          (TIME,UPAR,NPAR,DISP,VELO,JFLAG,IFLAG,

     &           RESULT,JDISP,JVELO)

C---- TO EXPORT * SUBROUTINE

      !DEC$ ATTRIBUTES DLLEXPORT,C::MATRIX_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     DISP   : Displacement vector between two force markers w.r.t. base marker. (Input)

C     VELO   : Velocity vector between two force markers w.r.t. base marker. (Input)

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 matrix force vector. (Output, Size : 6)

C     JDISP  : Returned force Jacobian of displacement. (Output, Size : 36)

C     JVELO  : Returned force Jacobian of velocity. (Output, Size : 36)

      DOUBLE PRECISION TIME, UPAR(*), DISP(*), VELO(*)

      INTEGER NPAR

      LOGICAL JFLAG, IFLAG

      DOUBLE PRECISION RESULT[REFERENCE](6)

      DOUBLE PRECISION JDISP[REFERENCE](36), JVELO[REFERENCE](36)

C---- USER STATEMENT

      DOUBLE PRECISION TK, TC, RK, RC

      INTEGER I

      TK = UPAR(1)

      TC = UPAR(2)

      RK = UPAR(3)

      RC = UPAR(4)

      IF ( JFLAG ) THEN

        DO I = 1, 36

          JDISP(I) = 0.0D0

          JVELO(I) = 0.0D0

        ENDDO

        JDISP(1) = -TK

        JDISP(8) = -TK

        JDISP(15) = -TK

        JDISP(22) = -RK

        JDISP(29) = -RK

        JDISP(36) = -RK

        JVELO(1) = -TC

        JVELO(8) = -TC

        JVELO(15) = -TC

        JVELO(22) = -RC

        JVELO(29) = -RC

        JVELO(36) = -RC

      ENDIF

      RESULT(1) = - TK * DISP(1) - TC * VELO(1)

      RESULT(2) = - TK * DISP(2) - TC * VELO(2)

      RESULT(3) = - TK * DISP(3) - TC * VELO(3)

      RESULT(4) = - RK * DISP(4) - RC * VELO(4)

      RESULT(5) = - RK * DISP(5) - RC * VELO(5)

      RESULT(6) = - RK * DISP(6) - RC * VELO(6)

      RETURN

      END