C---- SUB. MODAL_FORCE

      SUBROUTINE MODAL_FORCE

     &          (ID,TIME,UPAR,NPAR,IFBODY,POS,VEL,ACC,

     &           NMODE,NNODE,NMODALLOAD,MODALLOADS,

     &           JFLAG,IFLAG,RESULT)

      implicit none

C---- TO EXPORT * SUBROUTINE

      !DEC$ ATTRIBUTES DLLEXPORT,C::MODAL_FORCE

C---- INCLUDE SYSTEM CALL

      INCLUDE 'SYSCAL.F'

C---- DEFINE VARIABLES

c     Parameter Information

c     ID         : Modal force sequential identification. (Input)

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

c     UPAR       : Parameters defined by user. (Input)

c     NPAR       : Number of user parameters. (Input)

c     IFBODY     : RFLEX body sequential ID. (Input)

c     POS        : Position [1~3] and Orientation matrix [4~12] w.r.t Ground.InertiaMarker. (Input)

c     VEL        : Velocity vector w.r.t. Ground.InertiaMarker. (Input)

c     ACC        : Acceleration vector w.r.t Ground.InertiaMarker. (Input)

c     NMODE      : Number of selected mode. (Input)

c     NNODE      : Number of node. (Input)

c     NMODALLOAD : Number of selected modal load cases. (Input)

c     MODALLOADS : Modal force vector. (Input, Size : (6 + nmode) x nModalLoad)

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 modal froce vector. (Output, Size : 6 + nmode)

      INTEGER ID, NPAR, IFBODY

      DOUBLE PRECISION TIME, UPAR(*), POS(12), VEL(6), ACC(6)

      INTEGER NMODE, NNODE, NMODALLOAD

      DOUBLE PRECISION MODALLOADS(6+NMODE, NMODALLOAD)

      LOGICAL JFLAG, IFLAG

      DOUBLE PRECISION RESULT[REFERENCE](6+NMODE)

C----USER STATEMENT

      INTEGER i,j,ierr

      INTEGER mid(9)

      DOUBLE PRECISION RHO,Area(3),Cd,darea

      DOUBLE PRECISION MVel(3),VM,Fdir(3)

      DOUBLE PRECISION Fscale(3),dtmp(3),UCF

      LOGICAL eflag

C     INITIALIZE

      RHO = upar(1)

      Cd = upar(2)

      Area(1) = upar(3)

      Area(2) = upar(4)

      Area(3) = upar(5)

      mid(1) = int(upar(6))

      mid(2) = int (upar(7))

      mid(3) = int (upar(8))

      mid(4) = int (upar(9))

      mid(5) = int (upar(10))

      mid(6) = int (upar(11))

      mid(7) = int (upar(12))

      mid(8) = int (upar(13))

      mid(9) = int (upar(14))

      CALL RD_UCF(UCF)

      DO i=1, 6+nmode

        result(i) = 0.0d0

      ENDDO

      DO i=1, 9

        CALL sysary('TVEL',mid(i),1,MVel,3,eflag)

        VM = sqrt(MVel(1)*MVel(1)+MVel(2)*MVel(2)+MVel(3)*MVel(3))

        IF(VM .LT. 1.0d-17) THEN

          Fdir(1) = 0.0

          Fdir(2) = 0.0

          Fdir(3) = 0.0

        ELSE

          Fdir(1) = -MVel(1)/VM

          Fdir(2) = -MVel(2)/VM

          Fdir(3) = -MVel(3)/VM

        ENDIF

        IF(i.le.3) THEN

          darea = Area(1)

        ELSE IF (i.le.6) THEN

          darea = Area(2)

        ELSE

          darea = Area(3)

        ENDIF

        dtmp(1)=(Cd*RHO*VM*VM*darea/2)*Fdir(1)/UCF

        dtmp(2)=(Cd*RHO*VM*VM*darea/2)*Fdir(2)/UCF

        dtmp(3)=(Cd*RHO*VM*VM*darea/2)*Fdir(3)/UCF

        CALL ats(pos(4),dtmp,Fscale)

        DO j=1, 6+nmode

          result(j)=result(j)+

     &              ModalLoads(j,(i-1)*3+1)*Fscale(1)+

     &              ModalLoads(j,(i-1)*3+2)*Fscale(2)+

     &              ModalLoads(j,(i-1)*3+3)*Fscale(3)

        ENDDO

      ENDDO

      RETURN

      END

      SUBROUTINE ats(a,s,sp)

        implicit none

        DOUBLE PRECISION a(9),s(3),sp(3)

        sp(1)=A(1)*s(1)+A(2)*s(2)+A(3)*s(3)

        sp(2)=A(4)*s(1)+A(5)*s(2)+A(6)*s(3)

        sp(3)=A(7)*s(1)+A(8)*s(2)+A(9)*s(3)

      RETURN

      END