!#######################################################################
!     SUBROUTINE   --- sub_calc_density ---                          ###
!#######################################################################
!                                                                    ###
!     Subroutine for Calculating the Density of SPH Particles        ###
!                                                                    ###
!        Calculated Variable                                         ###
!           d(i) : Density of SPH particle                           ###
!                                                                    ###
!                                                         2009/07/29 ###
!                                                         2009/08/27 ###
!                                                         2011/06/27 ###
!#######################################################################

      subroutine sub_calc_density
!$ use omp_lib
      use mod_parameter
      use mod_variable

      implicit none

!==== Variables Used in This Subroutine ================================
      integer i
      integer j,jj
      real*8  wij
      real*8  w,dwdh,dphidh
      real*8  rij
      real*8  xi(3)
      real*8  h2i
      integer nbi_tot
      integer(4) :: trd
      real*8  omega_tmp
      real*8  zeta_tmp
!=======================================================================

      trd=0


!==== Calculating Density ==============================================
!$ call omp_set_num_threads(N_TRD)
!$OMP parallel default(none) &
!$OMP          private(i,h2i,xi,nbi_tot, &
!$OMP                  jj,j,wij,trd,rij,omega_tmp,zeta_tmp) &
!$OMP          shared(np,d,h,x,m,nbi_list,dis2_nbi_list,omega,zeta)

!$ trd=omp_get_thread_num()

!$OMP do
      do i=1,np

!        ---- Initialization ----
         d(i)=0.0d0
!        ------------------------

!        ---- Get the Neighbor List of i-th Particle -------------------
         h2i=2.0d0*h(i)
         xi(1)=x(1,i)
         xi(2)=x(2,i)
         xi(3)=x(3,i)
!##### 改良点：2*hi内だけの粒子をゲットすればよい ############
         call sub_get_neighbor_list(trd,xi,h2i,nbi_tot)
!#############################################################
!        ---------------------------------------------------------------

         do jj=1,nbi_tot
            j=nbi_list(jj,trd)
            wij=w(dis2_nbi_list(jj,trd),h(i))
            d(i)=d(i)+m(j)*wij
         end do

         !==== Calculation of Omega =================================
         omega_tmp=0.0d0
         do jj=1,nbi_tot
            j=nbi_list(jj,trd)
            rij=dsqrt(dis2_nbi_list(jj,trd))
            omega_tmp=omega_tmp+m(j)*dwdh(rij,h(i))
         end do
         omega(i)=1.0d0+h(i)/3.0d0/d(i)*omega_tmp
         !===========================================================

         !==== Calculation of Zeta ==================================
         zeta_tmp=0.0d0
         do jj=1,nbi_tot
            j=nbi_list(jj,trd)
            rij=dsqrt(dis2_nbi_list(jj,trd))
            zeta_tmp=zeta_tmp+m(j)*dphidh(rij,h(i))
         end do
         zeta(i)=-h(i)/3.0d0/d(i)*zeta_tmp
         !===========================================================

      end do
!$OMP end do
!$OMP end parallel

      return
      end