!#######################################################################
!     SUBROUTINE   --- sub_calc_smoothing_length ---                 ###
!#######################################################################
!                                                                    ###
!     Subroutine for Calculating Smoothing Length of SPH particle    ###
!                                                                    ###
!         Calculated Variables                                       ###
!             h(i)    ! Somoothing Length                            ###
!             nb(NPM) ! Number of Neighbor Particles                 ###
!                                                                    ###
!                                                         2009/07/24 ###
!                                                         2009/08/27 ###
!                                                         2011/06/23 ###
!#######################################################################

      subroutine sub_calc_smoothing_length
!$ use omp_lib
      use mod_parameter
      use mod_variable

      implicit none

!==== Variables Used in This Subroutine ================================
      integer i                     ! Particle Number
      real*8  c_nb,c_nb_damp        ! Dummy Variables
      real*8  xi(3)                 ! Position (x,y,z)
      real*8  h2i                   ! 2*h(i)
      integer nbi_tot               ! Number of Neighbor Particles
                                    !               for i-th Particle
      integer hi_ite                ! Number of Iteration
                                    !            for i-th Particle
      integer h_ite_tot             ! Total Number of Iteration
      real*8  hi_tmp_max            ! Dummy for Max h
      real*8  hi_tmp_min            ! Dummy for Min h
      integer(4) :: trd
!-----------------------------------------------------------------------

!==== Hernquist & Katz 1989 ApJS Eq.(2.17) =============================
      h_ite_tot=0
      h_ite_max=0
      trd=0

!$ call omp_set_num_threads(N_TRD)
!$OMP parallel default(none) &
!$OMP          private(i,hi_ite,hi_tmp_max,hi_tmp_min,c_nb_damp,c_nb, &
!$OMP                  h2i,xi,nbi_tot,trd) &
!$OMP          shared(np,x,nb,h,h_max,h_ite_tot,h_ite_max,nbi_list,dis2_nbi_list)

!$ trd=omp_get_thread_num()

!$OMP do reduction(+ : h_ite_tot) reduction(max : h_ite_max)
      do i=1,np

!       ---- Initialization ----
        hi_ite=1
        hi_tmp_max=1.0d30
        hi_tmp_min=0.0d0
!       ------------------------

!       ---- Start Iteration ----
  999   continue
!       ----------------------

!       ---- Original Version (Hernquist & Katz 1989) ----
!        c_nb=1.0d0+(dble(NB_STD)/dble(nb(i)))**(1.0d0/3.0d0)
!        h(i)=h(i)*0.5d0*c_nb
!       --------------------------------------------------

!       ---- Modified Version ----------------------------
        c_nb_damp=1.0d0+0.1d0*dble(hi_ite-1)
        c_nb=c_nb_damp+(dble(NB_STD)/dble(nb(i)))**(1.0d0/3.0d0)
        h(i)=h(i)*c_nb/(c_nb_damp+1.0d0)
!       --------------------------------------------------

        if(h(i) .lt. hi_tmp_min .or. h(i) .gt. hi_tmp_max) then
          h(i)=0.5d0*(hi_tmp_min+hi_tmp_max)
        endif

!==== Calculate the Number of Particles Within 2*h(i) ==================
        h2i=2.0d0*h(i)
        xi(1)=x(1,i)
        xi(2)=x(2,i)
        xi(3)=x(3,i)
        call sub_get_neighbor_number(trd,xi,h2i,nbi_tot)

!==== Judgement ========================================================
        nb(i)=nbi_tot

        if(nbi_tot .lt. NB_MIN) hi_tmp_min=h(i)
        if(nbi_tot .gt. NB_MAX) hi_tmp_max=h(i)
        if(FLAG_H_MAX == 1 .and. hi_tmp_min >= h_max) then
          h(i)=h_max
          go to 555
        endif

        if(hi_ite >= N_MAX_H_ITE) go to 555

        if(nbi_tot .gt. NB_MAX .or. nbi_tot .lt. NB_MIN) then
          hi_ite=hi_ite+1
          go to 999
        endif

  555   continue

!       ---- Cumlative and Max Number of Iteration ----
        h_ite_tot=h_ite_tot+hi_ite
        if(hi_ite .gt. h_ite_max) h_ite_max=hi_ite
!       -----------------------------------------------

      enddo
!$OMP end do
!$OMP end parallel

!==== Average Number of Iteration ======================================
      h_ite_ave=dble(h_ite_tot)/dble(np)


      return
      end