!#######################################################################
!     SUBROUTINE   --- sub_calc_COM ---                              ###
!#######################################################################
!                                                                    ###
!     Subroutine for Calculating Center of Mass                      ###
!                                                                    ###
!        Output Arguments                                            ###
!           xc(3) : Position of Center of Mass (x,y,z)               ###
!           vc(3) : Velocity of Center of Mass (x,y,z)               ###
!                                                                    ###
!                                                         2009/07/29 ###
!                                                         2009/08/27 ###
!                                                         2010/09/03 ###
!                                                         2011/06/23 ###
!#######################################################################

      subroutine sub_calc_COM(xc,vc)

      use mod_parameter
      use mod_variable

      implicit none

!==== Arguments ========================================================
!---- Output Arguments -------------------------------------------------
      real*8  xc(3)
      real*8  vc(3)

!==== Variables Used in This Subroutine ================================
      integer i           ! Particle Number
      integer k           ! Components (1-->x, 2-->y, 3-->z)
      real*8  m_tot       ! Total Mass
      real*8  mx_tot(3)   ! Summation of m(i)*x(k,i)
      real*8  mv_tot(3)   ! Summation of m(i)*v(k,i)
!=======================================================================

!     ---- Initialization ----
      m_tot=0.0d0
      mx_tot(1)=0.0d0
      mx_tot(2)=0.0d0
      mx_tot(3)=0.0d0
      mv_tot(1)=0.0d0
      mv_tot(2)=0.0d0
      mv_tot(3)=0.0d0
!     ------------------------

      do i=1,np
         m_tot=m_tot+m(i)
      enddo

      do i=1,np
         do k=1,3
            mx_tot(k)=mx_tot(k)+m(i)*x(k,i)
            mv_tot(k)=mv_tot(k)+m(i)*v(k,i)
         enddo
      enddo

      do k=1,3
         xc(k)=mx_tot(k)/m_tot
         vc(k)=mv_tot(k)/m_tot
      enddo


      return
      end