C---- SUB. NODAL_FORCE_EXT

      SUBROUTINE NODAL_FORCE_EXT

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

     &           JFLAG,IFLAG,RESULT)

C---- TO EXPORT * SUBROUTINE

      !DEC$ ATTRIBUTES DLLEXPORT,C::NODAL_FORCE_EXT

C---- INCLUDE SYSTEM CALL

      INCLUDE 'SYSCAL.F'

C---- DEFINE VARIABLES

C     Parameter Information

C     ID     : Nodal 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 : FFLEX Body sequential ID. (Input)

C     NODARR : Node sequential 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     IFLAG  : 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

c     local variable

      integer i,j,nodeid,errflg

      double precision Area, Cd, rho

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

      integer MKID(2), nMK

      Cd = upar(1)

      rho = upar(2)

      scale = upar(3)

      call rd_ucf(ucf)

      MKID(1) = 0

      MKID(2) = 0

      nMK = 0

      do i=1, nonde

         call get_fflex_nodetvel(ifbody,nodarr(i),MKID,nMK,vn,errflg)

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

         call get_fflex_nodeid(ifbody,nodarr(i),nodeid,errflg)

         if( nodeid .eq. 10000 .or. nodeid .eq. 10001) then

            Area = 100

         else if (nodeid .ge. 20000 .and. nodeid .lt. 30000) then

            Area = 200

         else Area = 400

         endif

         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

            endif

         enddo

      enddo

      RETURN

      END