!#######################################################################
!     SUBROUTINE   --- sub_calc_pot_energy_tree ---                  ###
!#######################################################################
!                                                                    ###
!     Subroutine for Calculating Gravitational Potential             ###
!     by Using Tree Method                                           ###
!                                                                    ###
!        Output Argument                                             ###
!           pot_eng : Total Gravitational Potential (Negative)       ###
!                                                                    ###
!                                                         2010/08/27 ###
!                                                         2011/06/23 ###
!#######################################################################

      subroutine sub_calc_pot_energy_tree(pot_eng)

!$ use omp_lib
      use mod_parameter
      use mod_variable

      implicit none

!==== Output Argument ==================================================
      real*8  pot_eng

!==== Variables Used in This Subroutine ================================
      integer i
      integer j
      integer k
      real*8  rij2
      real*8  rij
      integer gt_tot
      real*8  mc_rij_rev
      real*8  xi(3)
      integer jj
      real*8  hij
      real*8  f_pot_s
      real*8  theta_g2
      integer(4) :: trd
!=======================================================================

      theta_g2=THETA_GRAV2_ENG
      pot_eng=0.0d0
      trd=0

!$ call omp_set_num_threads(N_TRD)
!$OMP parallel default(none) &
!$OMP          private(i,k,xi,gt_tot,j,rij2,rij,jj, &
!$OMP                  hij,mc_rij_rev,trd) &
!$OMP          shared(np,x,h,hg,pot_eng,m,theta_g2,id_nb,mc_nb,xc_nb,dis2_nb)

!$ trd=omp_get_thread_num()

!$OMP do reduction(+ : pot_eng)
      do i=1,np
         do k=1,3
            xi(k)=x(k,i)
         enddo

         call sub_get_grav_tree_list(trd,theta_g2,xi,h(i),gt_tot)

         do j=1,gt_tot
            rij2=dis2_nb(j,trd)
            rij=dsqrt(rij2)
            if(id_nb(j,trd) .ne. -1) then
               jj=id_nb(j,trd)
               mc_rij_rev=mc_nb(j,trd)*0.5d0*(f_pot_s(rij,h(i))/h(i)+f_pot_s(rij,h(jj))/h(jj))
            else
               mc_rij_rev=mc_nb(j,trd)/rij
            endif
            pot_eng=pot_eng-0.5d0*GC*m(i)*mc_rij_rev
         enddo
      enddo
!$OMP end do
!$OMP end parallel


      return
      end


!#######################################################################
!     SUBROUTINE   --- sub_calc_pot_energy_direct ---                ###
!#######################################################################
!                                                                    ###
!     Subroutine for Calculating Gravitational Potential             ###
!     by Using Direct Method                                         ###
!                                                                    ###
!        Output Argument                                             ###
!           pot_eng : Total Gravitational Potential (Negative)       ###
!                                                                    ###
!                                                         2010/09/03 ###
!                                                         2011/06/23 ###
!#######################################################################

subroutine sub_calc_pot_energy_direct(pot_eng)

!$ use omp_lib
  use mod_parameter
  use mod_variable

  implicit none

!==== Output Argument ==================================================
  real(8), intent(OUT) :: pot_eng

!==== Variables Used in This Subroutine ================================
  integer(4) :: i
  integer(4) :: j
  real(8)    :: rij2
  real(8)    :: rij
  real(8)    :: h2_max
  real(8)    :: m_rij_rev
  real(8)    :: f_pot_s
  integer(4) :: trd
!=======================================================================

  pot_eng=0.0d0

!$ call omp_set_num_threads(N_TRD)
!$OMP parallel default(none) &
!$OMP          private(i,j,rij2,rij,h2_max,m_rij_rev,trd) &
!$OMP          shared(np,x,h,m,pot_eng)

!$ trd=omp_get_thread_num()

!$OMP do reduction(+ : pot_eng)
  do i=1,np
     do j=1,np
        rij2=(x(1,i)-x(1,j))**2 &
            +(x(2,i)-x(2,j))**2 &
            +(x(3,i)-x(3,j))**2
        rij=dsqrt(rij2)
        h2_max=dmax1(2.0d0*h(i),2.0d0*h(j))
        if(rij >= h2_max) then
           m_rij_rev=m(j)/rij
        else
           m_rij_rev=m(j)*0.5d0*( f_pot_s(rij,h(i))/h(i) &
                               & +f_pot_s(rij,h(j))/h(j) )
        endif
        pot_eng=pot_eng-0.5d0*GC*m(i)*m_rij_rev
     end do
  end do
!$OMP end do
!$OMP end parallel

  return
end subroutine sub_calc_pot_energy_direct


!#######################################################################
!     FUNCTION   --- f_pot_s ---                                     ###
!#######################################################################
!                                                                    ###
!     Function for Calculating Gravitational Potential               ###
!     between two SPH Particles                                      ###
!                                                                    ###
!        Input Arguments                                             ###
!           r  : Distance between two SPH Particles                  ###
!           h  : Average Smoothing Length between two SPH Particles  ###
!                                                                    ###
!                                                         2011/06/23 ###
!#######################################################################

      real*8 function f_pot_s(r,h)
      implicit none

!==== Input Arguments ==================================================
      real*8  r,h
!==== Variables Used in This Function ==================================
      real*8  s,s2,s3,s4,s5
!=======================================================================

      s=r/h

      if (s .ge. 2.0d0) then
         f_pot_s=1.0d0/s
      else if (s .ge. 1.0d0) then
            s2=s**2
            s3=s2*s
            s4=s2*s2
            s5=s3*s2
         f_pot_s=-1.0d0/15.0d0/s+8.0d0/5.0d0-4.0d0/3.0d0*s2+s3 &
                 -3.0d0/10.0d0*s4+1.0d0/30.0d0*s5
      else
            s2=s**2
            s4=s2*s2
            s5=s4*s
         f_pot_s=7.0d0/5.0d0-2.0d0/3.0d0*s2+3.0d0/10.0d0*s4 &
                 -1.0d0/10.0d0*s5
      end if


      return
      end