!#######################################################################
!     SUBROUTINE   --- sub_calc_time_step ---                        ###
!#######################################################################
!                                                                    ###
!     Subroutine for Determining Time Step                           ###
!     Determination of time step is based on Monaghan (1989).        ###
!                                                                    ###
!        Input Arguments                                             ###
!           nu_max : Max Artificial Viscosity                        ###
!        Output Argument                                             ###
!           dt     : Time Step                                       ###
!                                                                    ###
!                                                         2009/07/29 ###
!                                                         2009/08/27 ###
!                                                         2011/06/26 ###
!#######################################################################

subroutine sub_calc_time_step(dt)

  use mod_parameter
  use mod_variable

  implicit none

!==== Arguments ========================================================
!---- Output Arguments -------------------------------------------------
  real(8), intent(OUT) :: dt

!==== Variables Used in This Subroutine ================================
  integer(4) :: i           ! Particle Number
  real(8)    :: dt_v_tmp,dt_h_tmp
  real(8)    :: dt_v,dt_h,dt_e
  real(8)    :: a2          ! Acceleration Squared
!=======================================================================

  real(8), parameter :: EDT=0.1d0

  dt_v=1.0d30
  dt_h=1.0d30
  dt_e=1.0d30

  do i=1,np

     !---- Limit by Viscosity ------------------------------------------
     dt_v_tmp=h(i)/(sv(i)+1.2d0*(ALP_VIS*sv(i)+BET_VIS*nu_max(i)))
     dt_v=dmin1(dt_v,dt_v_tmp)

     !---- Limit by Smoothing Length -----------------------------------
     a2=a(1,i)**2+a(2,i)**2+a(3,i)**2
     if (dabs(a2) > 1.0d-60) then
        dt_h_tmp=dsqrt(h(i)/dsqrt(a2))
        dt_h=dmin1(dt_h,dt_h_tmp)
     end if

  end do

  dt=dmin1(CDT*dt_v,CDT*dt_h,EDT*dt_e)

  return
end subroutine sub_calc_time_step