!#######################################################################
!     SUBROUTINE   --- sub_next_step ---                             ###
!#######################################################################
!                                                                    ###
!     Subroutine for Calculating Positions, Velocities, and          ###
!     Specific Internal Energy of All SPH Particles at the Next Step ###
!                                                                    ###
!        Input Arguments                                             ###
!           n_solve     ! Iteration Count                            ###
!           dt          ! Time Step                                  ###
!                                                                    ###
!        Calculated Variables (Common Blocks)                        ###
!           x(3,NPM)   : Position (x,y,z)                            ###
!           v(3,NPM)   : Velocity (x,y,z)                            ###
!           e(NPM)     : Specific Internal Energy                    ###
!                                                                    ###
!                                                         2009/09/02 ###
!                                                         2011/06/23 ###
!                                                         2011/09/21 ###
!#######################################################################

subroutine sub_next_step(n_step,n_solve,dt)

  use mod_parameter
  use mod_variable

  implicit none

!==== Input Arguments ==================================================
  integer(4), intent(IN) :: n_step
  integer(4), intent(IN) :: n_solve
  real(8),    intent(IN) :: dt

!==== Variables Used in This Subroutine ================================
  integer(4) :: i             ! Particle Number
  integer(4) :: k             ! Component (x,y,z)
  real(8)    :: dt2_half      ! 0.5*dt**2
  real(8)    :: de
  real(8)    :: dedt2
  real(8)    :: v2, v1, v1_max
!=======================================================================

  dt2_half=0.5d0*dt*dt

!==== Case for 1st-Order Euler Method ==================================
  if (FLAG_SOLVE .eq. 0) then
     do i=1,np
        do k=1,3
           x(k,i)=x(k,i)+v(k,i)*dt+a(k,i)*dt2_half
           v(k,i)=v(k,i)+a(k,i)*dt
        enddo
        de=dedt(i)*dt
        if(dedt(i) < 0.0d0 .and. 2.0d0*dabs(de) > e(i)) then
           if(FLAG_ALERT == 1) call sub_alert_negative_e(n_step,i,n_solve,e(i))
           e(i)=0.5d0*e(i)
        else
           e(i)=e(i)+de
        endif
     enddo
!    %%%%%%%%%% Relaxation Limit %%%%%%%%%%
     if (FLAG_RELAX .eq. 1) then
        do i=1,np
           v2=v(1,i)**2+v(2,i)**2+v(3,i)**2
           v1=sqrt(v2)
           v1_max=max(CSV_RELAX*sv(i),VMAX_RELAX)
           if(v1 > v1_max) then
              do k=1,3
                 v(k,i)=v(k,i)/v1*v1_max
              enddo
           endif
        enddo
     endif
!    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  endif

!==== Case for 2nd-Order PEC Method ====================================
  if (FLAG_SOLVE .eq. 1) then
     if (n_solve .eq. 1) then
        do i=1,np
           do k=1,3
              x(k,i)=x(k,i)+v(k,i)*dt+a(k,i)*dt2_half
              v(k,i)=v(k,i)+a(k,i)*dt
              a1(k,i)=a(k,i)
           enddo
           de=dedt(i)*dt
           if(dedt(i) < 0.0d0 .and. 2.0d0*dabs(de) > e(i)) then
              if(FLAG_ALERT == 1) call sub_alert_negative_e(n_step,i,n_solve,e(i))
              e(i)=0.5d0*e(i)
              dedt1(i)=e(i)/dt
           else
              e(i)=e(i)+de
              dedt1(i)=dedt(i)
           endif
        enddo
     else
        do i=1,np
           do k=1,3
!fastfastfastfastfastfast
              x(k,i)=x(k,i)+PEC_A*(a(k,i)-a1(k,i))*dt2_half
              v(k,i)=v(k,i)+PEC_B*(a(k,i)-a1(k,i))*dt
!fastfastfastfastfastfast
           end do
           dedt2=dedt(i)-dedt1(i)
           de=PEC_C*dedt2*dt
           if(dedt2 < 0.0d0 .and. 2.0d0*dabs(de) > e(i)) then
              if(FLAG_ALERT == 1) call sub_alert_negative_e(n_step,i,n_solve,e(i))
              e(i)=0.5d0*e(i)
           else
              e(i)=e(i)+de
           endif
        end do
     endif
  endif


  return
end subroutine sub_next_step


!#######################################################################
!     SUBROUTINE   --- sub_next_step_e_table ---                     ###
!#######################################################################
!                                                                    ###
!     Subroutine for Calculating Positions, Velocities, and          ###
!     Specific Internal Energy of All SPH Particles at the Next Step ###
!                                                                    ###
!        Input Arguments                                             ###
!           n_solve     ! Iteration Count                            ###
!           dt          ! Time Step                                  ###
!                                                                    ###
!        Calculated Variables (Common Blocks)                        ###
!           x(3,NPM)   : Position (x,y,z)                            ###
!           v(3,NPM)   : Velocity (x,y,z)                            ###
!           e(NPM)     : Specific Internal Energy                    ###
!                                                                    ###
!                                                         2009/09/02 ###
!                                                         2011/06/23 ###
!                                                         2011/09/21 ###
!#######################################################################

subroutine sub_next_step_e_table(n_solve,dt, &
&                        i_tb_max,dr_tb,r_tb,e_tb)

  use mod_parameter
  use mod_variable

  implicit none

!==== Input Arguments ==================================================
  integer(4), intent(IN) :: n_solve
  real(8),    intent(IN) :: dt

!==== Variables Used in This Subroutine ================================
  integer(4) :: i             ! Particle Number
  integer(4) :: k             ! Component (x,y,z)
  real(8)    :: dt2_half      ! 0.5*dt**2
!=======================================================================

  real(8) :: dr_tb
  integer(4) :: i_tb_max
  real(8) ::  r_tb(0:1000),e_tb(0:1000)
  real(8) :: r,r2
  integer(4) :: i_tb

  dt2_half=0.5d0*dt*dt

!==== Case for 1st-Order Euler Method ==================================
  if (FLAG_SOLVE .eq. 0) then
     do i=1,np
        do k=1,3
           x(k,i)=x(k,i)+v(k,i)*dt+a(k,i)*dt2_half
           v(k,i)=v(k,i)+a(k,i)*dt
        enddo
        !---- テーブルを使って内部エネルギーを与える ----
        r2=x(1,i)**2+x(2,i)**2+x(3,i)**2
        r=dsqrt(r2)
        i_tb=int(r/dr_tb)
        if (i_tb .ge. i_tb_max) then
           e(i)=0.0d0
        else
           e(i)=e_tb(i_tb) &
           &             +(e_tb(i_tb+1)-e_tb(i_tb))/dr_tb*(r-r_tb(i_tb))
        endif
        !------------------------------------------------
!        de=dedt(i)*dt
!        if(dedt(i) < 0.0d0 .and. 2.0d0*dabs(de) > e(i)) then
!           call sub_alert_negative_e(n_step,i,n_solve,e(i))
!           e(i)=0.5d0*e(i)
!        else
!           e(i)=e(i)+de
!        endif
     enddo
  endif

!==== Case for 2nd-Order PEC Method ====================================
  if (FLAG_SOLVE .eq. 1) then
     if (n_solve .eq. 1) then
        do i=1,np
           do k=1,3
              x(k,i)=x(k,i)+v(k,i)*dt+a(k,i)*dt2_half
              v(k,i)=v(k,i)+a(k,i)*dt
              a1(k,i)=a(k,i)
           enddo
        !---- テーブルを使って内部エネルギーを与える ----
        r2=x(1,i)**2+x(2,i)**2+x(3,i)**2
        r=dsqrt(r2)
        i_tb=int(r/dr_tb)
        if (i_tb .ge. i_tb_max) then
           e(i)=0.0d0
        else
           e(i)=e_tb(i_tb) &
           &             +(e_tb(i_tb+1)-e_tb(i_tb))/dr_tb*(r-r_tb(i_tb))
        endif
        !------------------------------------------------
!           de=dedt(i)*dt
!           if(dedt(i) < 0.0d0 .and. 2.0d0*dabs(de) > e(i)) then
!              call sub_alert_negative_e(n_step,i,n_solve,e(i))
!              e(i)=0.5d0*e(i)
!              dedt1(i)=e(i)/dt
!           else
!              e(i)=e(i)+de
              dedt1(i)=dedt(i)
!           endif
        enddo
     else
        do i=1,np
           do k=1,3
!fastfastfastfastfastfast
              x(k,i)=x(k,i)+PEC_A*(a(k,i)-a1(k,i))*dt2_half
              v(k,i)=v(k,i)+PEC_B*(a(k,i)-a1(k,i))*dt
!fastfastfastfastfastfast
           end do
        !---- テーブルを使って内部エネルギーを与える ----
        r2=x(1,i)**2+x(2,i)**2+x(3,i)**2
        r=dsqrt(r2)
        i_tb=int(r/dr_tb)
        if (i_tb .ge. i_tb_max) then
           e(i)=0.0d0
        else
           e(i)=e_tb(i_tb) &
           &             +(e_tb(i_tb+1)-e_tb(i_tb))/dr_tb*(r-r_tb(i_tb))
        endif
        !------------------------------------------------
!           dedt2=dedt(i)-dedt1(i)
!           de=PEC_C*dedt2*dt
!           if(dedt2 < 0.0d0 .and. 2.0d0*dabs(de) > e(i)) then
!              call sub_alert_negative_e(n_step,i,n_solve,e(i))
!              e(i)=0.5d0*e(i)
!           else
!              e(i)=e(i)+de
!           endif
        end do
     endif
  endif


  return
end subroutine sub_next_step_e_table


!#######################################################################
!     SUBROUTINE   --- sub_alert_negative_e ---                      ###
!#######################################################################
!                                                                    ###
!     Subroutine for Output of Alert                                 ###
!                                                                    ###
!                                                         2011/09/21 ###
!#######################################################################

subroutine sub_alert_negative_e(n_step,i,n_solve,e)
  implicit none

!==== Input Arguments ==================================================
  integer(4), intent(IN) :: n_step
  integer(4), intent(IN) :: i
  integer(4), intent(IN) :: n_solve
  real(8),    intent(IN) :: e

!==== Output Alert =====================================================
  write(92,*) '===== Negative Internal Energy ====='
  write(92,'(a7,$)')    '   STEP'
  write(92,'(a7,$)')    '   ID #'
  write(92,'(a7,$)')    '    PEC'
  write(92,'(a12,$)')   '    OLD e(i)'
  write(92,'(a12)')     '    NEW e(i)'
  write(92,'(i7,$)')    n_step
  write(92,'(i7,$)')    i
  write(92,'(i7,$)')    n_solve
  write(92,'(e12.4,$)') e
  write(92,'(e12.4)')   0.5d0*e
  write(92,*) ' '
  call flush(92)


  return
end subroutine sub_alert_negative_e