C---- SUB. SCREW_FORCE

      SUBROUTINE SCREW_FORCE

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

C---- TO EXPORT * SUBROUTINE

      !DEC$ ATTRIBUTES DLLEXPORT,C::SCREW_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     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 screw force vector. (Output, Size : 6)

      DOUBLE PRECISION TIME, UPAR(*)

      INTEGER NPAR

      LOGICAL JFLAG, IFLAG

      DOUBLE PRECISION RESULT[REFERENCE](6)

C---- USER STATEMENT

C---- USER DECLARED VARIABLES

      INTEGER MKID(2)

      DOUBLE PRECISION VALUE1, VALUE2, FORCE, ISIGN, L, MU, K, C

      LOGICAL ERRFLG

      L = 200.0D0

      MU = 0.3D0

      K = 1000.0D0

      C = 10.0D0

C--- COMPUTE A NORMAL FORCE

C---- AFTER CASTING DOUBLE PRECISION TYPE ‘UPAR’INTO INTEGER TYPE ‘MKID’

C---- ASSIGNS ‘UPAR(1), (3)’INTO 'MKID(1),(2)'.

      MKID(1) = INT(UPAR(1))

      MKID(2) = INT(UPAR(3))

C---- BY CALLING ‘SYSFNC’, USER CAN OBTAIN 'VALUE1' WHICH MEANS 'DY'.

      CALL SYSFNC('DY',MKID,2,VALUE1,ERRFLG)

C---- LIKEWISE, USER CAN OBTAIN 'VALUE2' WHICH MEANS 'VY'.

      CALL SYSFNC('VY',MKID,2,VALUE2,ERRFLG)

C---- WHEN THE WHEEL PENETRATES INTO THE GROUND, THE ‘RESULT(2)’WHICH MEANS ‘FY’HAS A VALUE.

      RESULT(2) = -K*(VALUE1-L) - C*VALUE2

      IF ( RESULT(2) .LE. 0.0D0 ) RESULT(2) = 0.0D0

C--- COMPUTE A FRICTION FORCE

      CALL SYSFNC('VX',MKID,2,VALUE1,ERRFLG)

      CALL SYSFNC('WZ',MKID,2,VALUE2,ERRFLG)

      IF ( VALUE1+VALUE2*L .LT. -0.01D0 ) THEN

        ISIGN = -1.0D0

      ELSE IF ( VALUE1+VALUE2*L .GT. 0.01D0 ) THEN

        ISIGN = 1.0D0

      ELSE

        ISIGN = 0.0D0

      ENDIF

      MKID(1) = INT(UPAR(1))

      MKID(2) = INT(UPAR(2))

      CALL SYSFNC('FY',MKID,2,FORCE,ERRFLG)

      RESULT(1) = -MU*DABS(FORCE)*ISIGN

C---- COMPUTE A TORQUE GENERATED BY THE FRICTION FORCE

      CALL SYSFNC('FX',MKID,2,FORCE,ERRFLG)

      RESULT(6) = L*DABS(FORCE)

      RETURN

      END