C---- SUB. MODAL_FORCE_EXT

      SUBROUTINE MODAL_FORCE_EXT

     &          (ID,TIME,UPAR,NPAR,IFBODY,NODARR,NONDE,

     &           JFLAG,IFLAG,RESULT)

C---- TO EXPORT * SUBROUTINE

      !DEC$ ATTRIBUTES DLLEXPORT,C::MODAL_FORCE_EXT

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     NODARR : Node ID array of input node set. (Input)

C     NONDE  : Number of node of node set. (Input)

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

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

C     RESULT : Returned nodal force vector. Acting point of the nodal force is each center of each node.

C              Reference frame of each force vector must be the ground inertia marker. (Output, Size: nonde * 6)

      DOUBLE PRECISION TIME, UPAR(*)

      INTEGER ID, NPAR, IFBODY, NODARR(*), NONDE

      LOGICAL JFLAG, IFLAG

      DOUBLE PRECISION RESULT(*)

C---- USER STATEMENT

      integer i,j,errflg

      double precision Area, Cd, rho

      double precision vn(3), vn_mag, ucf, scale

      integer MKID(2), nMK

      integer nodeseq

      Cd = upar(2)

      rho = upar(3)

      scale = upar(4)

      call rd_ucf(ucf)

      MKID(1) = 0

      MKID(2) = 0

      nMK = 0

      if(nodarr(i) .eq. 50004 .or. nodarr(i) .eq. 50007) then

        Area = 100

      else if (nodarr(i) .lt. 50010 .and. nodarr(i) .gt. 50075 )  then

        Area = 200

      else if(nodarr(i) .gt. 50076)  then

        Area = 400

      endif 

      do i=1, nonde

        call get_rflex_nodeseqid(ifbody,nodarr(i),nodeseq,errflg)

        call get_rflex_nodetvel(ifbody,nodeseq,MKID,nMK,vn,errflg)

        vn_mag = dsqrt(vn(1)**2+vn(2)**2+vn(3)**2)

        do j=1, 3

          if(vn_mag .eq. 0.0d0) then

            result((i-1)*6+j) = 0.0d0

          else

            result((i-1)*6+j)=-1.0d0/2.0d0*Cd*rho*vn_mag*vn_mag*Area

     &                         *(vn(j)/vn_mag) ! direction

     &                         /ucf            ! making force unit

     &                         *scale          ! an amplifier

            result((i-1)*6+j+3) = 0.0;

          endif

        enddo

      enddo

      RETURN

      END