!#######################################################################
!     SUBROUTINE   --- sub_calc_Angular_Momentum ---                 ###
!#######################################################################
!                                                                    ###
!     Subroutine for Calculating Total Angular Momentum              ###
!                                                                    ###
!        Input Arguments                                             ###
!           xc(3)   : Position of Center of Mass (x,y,z)             ###
!           vc(3)   : Velocity of Center of Mass (x,y,z)             ###
!        Output Argument                                             ###
!           amom(3) : Total Angular Momentum (x,y,z)                 ###
!                                                                    ###
!                                                         2010/09/22 ###
!                                                         2011/06/21 ###
!#######################################################################

      subroutine sub_calc_angular_momentum(xc,vc,amom)

      use mod_parameter
      use mod_variable

      implicit none

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

!==== Variables Used in This Subroutine ================================
      integer i           ! Particle Number
      integer k           ! Components (1-->x, 2-->y, 3-->z)
      real*8  xic(3)      ! x(k,i)-xc(k)
      real*8  vic(3)      ! v(k,i)-vc(k)
      real*8  damom(3)    ! Angular Momentum of Each SPH Particle
!=======================================================================

!     ---- Initialization ----
      amom(1)=0.0d0
      amom(2)=0.0d0
      amom(3)=0.0d0
!     ------------------------

      do i=1,np
         do k=1,3
            xic(k)=x(k,i)-xc(k)
            vic(k)=v(k,i)-vc(k)
         enddo
         damom(1)=m(i)*(vic(2)*xic(3)-vic(3)*xic(2))
         damom(2)=m(i)*(vic(3)*xic(1)-vic(1)*xic(3))
         damom(3)=m(i)*(vic(1)*xic(2)-vic(2)*xic(1))
         do k=1,3
            amom(k)=amom(k)+damom(k)
         enddo
      enddo


      return
      end