!#######################################################################
!     SUBROUTINE   --- sub_calc_accel ---                            ###
!#######################################################################
!                                                                    ###
!     Subroutine for Calculating the Acceleration and de/dt          ###
!                                                                    ###
!        Output Arguments                                            ###
!           nu_max  : Max Viscosity Used for Time Step Determination ###
!                                                                    ###
!                                                         2009/07/29 ###
!                                                         2009/08/27 ###
!                                                         2011/06/27 ###
!#######################################################################

subroutine sub_calc_accel
!$ use omp_lib
  use mod_parameter
  use mod_variable

  implicit none

!==== Variables Used in This Subroutine ================================
  integer(4) ::i
  real(8) :: d2i
  integer(4) ::j,jj
  real(8) :: xi(3)
  real(8) :: h2i
  integer(4) ::nbi_tot
  real(8) :: dx(3)
  real(8) :: dv(3)
  real(8) :: dis2
  real(8) :: dw
  real(8) :: vx
  real(8) :: hij
  real(8) :: nu,cs,q_half
  real(8) :: pddij
  real(8) :: mpddij
  real(8) :: dedt_temp
  integer(4) :: trd
  real(8) :: tmp,tmp_i,tmp_j
  real(8) :: dwi,dwj
  real(8) :: pdd_omega_i,pdd_omega_j
!-----------------------------------------------------------------------

  nu_max=0.0d0
  trd=0

!$ call omp_set_num_threads(N_TRD)
!$OMP parallel default(none) &
!$OMP          private(i,h2i,xi,nbi_tot, &
!$OMP                  jj,j,dx,dv,dis2,vx,hij,nu,cs,q_half, &
!$OMP                  pddij,mpddij,dedt_temp,trd,d2i,tmp,tmp_i,tmp_j,dwi,dwj,pdd_omega_i,pdd_omega_j) &
!$OMP          shared(np,h,x,v,nu_max,sv,d,pdd,p,m,e,dedt,a,nbi_list,dis2_nbi_list,omega)

!$ trd=omp_get_thread_num()

!---- Initialization ---------------------------------------------------
!$OMP do
  do i=1,np
     a(1:3,i)=0.0d0
     dedt(i)=0.0d0
     nu_max(i)=0.0d0
  enddo
!$OMP end do

!---- Pre Calculation --------------------------------------------------
!$OMP do
  do i=1,np
     d2i=d(i)**2
     pdd(i)=p(i)/d2i
  enddo
!$OMP end do

!==== Calculating Acceleration and de/dt ===============================
!$OMP do
  do i=1,np
     !---- Get the Neighbor List of i-th Particle ----------------------
     h2i=2.0d0*h(i)
     xi(1:3)=x(1:3,i)
     call sub_get_neighbor_list(trd,xi,h2i,nbi_tot)

     do jj=1,nbi_tot

        j=nbi_list(jj,trd)

        if (j .eq. i) go to 100

        dx(1:3)=x(1:3,i)-x(1:3,j)
        dv(1:3)=v(1:3,i)-v(1:3,j)

        dis2=dx(1)**2+dx(2)**2+dx(3)**2

        !---- Artificial Viscosity ----
        vx=dv(1)*dx(1)+dv(2)*dx(2)+dv(3)*dx(3)
        if ( vx .lt. 0.0d0 ) then
           hij=(h(i)+h(j))*0.5d0
           nu=vx/hij/(dis2/hij/hij+ETA_VIS2)
           nu_max(i)=dmax1(nu_max(i),dabs(nu))
           cs=(sv(i)+sv(j))/2.0d0
           q_half=(-ALP_VIS*nu*cs+BET_VIS*nu*nu)/(d(i)+d(j))
        else
           q_half=0.0d0
        end if
        !---- Relaxation ----
!        if (FLAG_RELAX == 1) q_half = 0.0e0_8
        !--------------------
        !------------------------------

        dwi=dw(dis2,h(i))
        dwj=dw(dis2,h(j))
        pdd_omega_i=pdd(i)/omega(i)
        pdd_omega_j=pdd(j)/omega(j)

        tmp_i = (pdd_omega_i+q_half)*dwi
        tmp_j = (pdd_omega_j+q_half)*dwj
        tmp = m(j)*(tmp_i+tmp_j)

        a(1:3,i)=a(1:3,i)-tmp*dx(1:3)
        dedt(i)=dedt(i) &
             & +m(j)*(pdd_omega_i*dwi+0.5d0*q_half*(dwi+dwj))*vx

100     continue

     enddo


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

  return
end subroutine sub_calc_accel