subroutine group(fm,lm,ifirst,ilast,
     &  x, polar, momen, rhoGrid, rho, dr, cw, ccw, Gforce)
c     ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      include "./linearad.i"
c     ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer fm,lm,ifirst,ilast, rhoGrid
      double precision x(fm:lm,1:4)
      double precision polar, momen
      double precision rho(1:rhoGrid,1:6)
      double precision dr
      integer cw, ccw
      double precision Gforce(fm:lm,1:2)

      integer ie, je, totaln
      double precision vector(1:2),scalar
      double precision cm(1:2), dis(1:2), torq, torq_ind
      double precision dist, max_dist, momr, momt, speed
c     ******************************************************************
c     Calculate Polarity
c     ******************************************************************
      totaln=(ilast-ifirst+1)
      vector(1)=0.d0
      vector(2)=0.d0
      scalar=0.d0
      do ie=ifirst, ilast
         vector(1)=vector(1)+x(ie,3)
         vector(2)=vector(2)+x(ie,4)
         scalar=scalar+( x(ie,3)**2+x(ie,4)**2 )**0.5d0
      enddo
      vector(1)=vector(1)/totaln
      vector(2)=vector(2)/totaln
      scalar=scalar/totaln
      polar=( vector(1)**2+vector(2)**2 )**0.5d0/scalar
c     ******************************************************************
c     Calculate Momentum
c     ******************************************************************
c     Center of Mass
      cw = 0
      ccw = 0
      cm(1)=0.d0
      cm(2)=0.d0
      do ie=ifirst,ilast
         cm(1)=cm(1)+x(ie,1)
         cm(2)=cm(2)+x(ie,2)
      enddo
      cm(1)=cm(1)/totaln
      cm(2)=cm(2)/totaln
c     Momentum
      torq=0.d0
      scalar=0.d0
      max_dist=0.d0
      do ie=ifirst,ilast
         dis(1)=x(ie,1)-cm(1)
         dis(2)=x(ie,2)-cm(2)
         dist = (dis(1)**2+dis(2)**2)**0.5d0
         if(dist>max_dist) then
            max_dist = dist
         endif
         torq_ind = ( dis(1)*x(ie,4)-dis(2)*x(ie,3) )
         if(torq_ind.gt.0.d0.and.dist.le.10.d0*la) then
            ccw = ccw+1
         else if(torq_ind.lt.0.d0.and.dist.le.10.d0*la) then
            cw = cw+1
         endif
         torq=torq+torq_ind
         scalar=scalar+ ( dis(1)**2+dis(2)**2 )**0.5d0
     &                 *( x(ie,3)**2+x(ie,4)**2 )**0.5d0
      enddo
      torq=torq/totaln
      scalar=scalar/totaln
      momen=torq/scalar

      speed = sqrt(alpha/beta)
      if(dr.le.0.0001) then
         dr = max_dist/dble(rhoGrid)*1.5d0
      endif
      do je=1,rhoGrid
         rho(je,1)=dble(je)*dr
         rho(je,2)=0.d0
         rho(je,3)=0.d0
         rho(je,4)=0.d0
         rho(je,5)=0.d0
         rho(je,6)=0.d0
      enddo
      do ie=ifirst,ilast
         dis(1)=x(ie,1)-cm(1)
         dis(2)=x(ie,2)-cm(2)
         dist = (dis(1)**2+dis(2)**2)**0.5d0
         momr=(x(ie,3)*dis(1)+x(ie,4)*dis(2))/dist
         momt=(x(ie,3)*dis(2)-x(ie,4)*dis(1))/dist
         je = int(dist/dr)+1
         if(je.le.rhoGrid) then
            rho(je,2)=rho(je,2)+1
            rho(je,3)=rho(je,3)+momr
            rho(je,4)=rho(je,4)+momt
            rho(je,5)=rho(je,5)+(abs(momt)-speed)**2
     &                         +(momr-0.d0)**2
            rho(je,6)=rho(je,6)+Gforce(ie,1)**2
     &                         +Gforce(ie,2)**2
         endif
      enddo

      return
      end

c*************************************************************************