!-----------------------------------------------------------------------------
!
!
                         SUBROUTINE dt_comp
!
!
!-----------------------------------------------------------------------------
!
!
!     Subroutine to compute the DT integration interval, FLY_2.0 Sept. 2001
!
!
!-----------------------------------------------------------------------------

	 USE fly_h 
	implicit none
  	 INCLUDE 'mpif.h'

!   Declaration of local variables.
!   -------------------------------
	REAL(KIND=8) :: dt_act,t_new,a_guess,z_new,max_dt_acc,mod_a
	REAL(KIND=8) ::max_acc_local,mod_acc,mod_vel,max_mod_acc,max_mod_vel,max_modg_acc,max_modg_vel
	REAL(KIND=8) :: ttn2,ap2,a_dot2,bp2
	INTEGER	     :: p,ibod,ibod2,ibod3
	REAL(KIND=8) :: dt_guess
!-----------------------------------------------------------------------------
!   max acceleration component in all the system
!-----------------------------------------------------------------------------
	
	max_acc_local=0.0
	mod_vel=0.
	max_mod_acc=0.
	max_mod_vel=0.
	max_modg_acc=0.
	max_modg_vel=0.
	
	IF(dt_var) THEN
        ttn2=tnow+dtime2
        ap2=f_ap/ttn2
        a_dot2 = om_sum+Omega_l*(ttn2**three_o_alpha)
        bp2 = 1./(alpha2*(hubble*hubble)*a_dot2*(ttn2**2.))
        DO p=1,nb_res_loc(me+1)
	  mod_acc=abs(acc(1,p)*acc(1,p) + &
               acc(2,p)*acc(2,p)+acc(3,p)*acc(3,p))
	  mod_vel=abs(vel(1,p)*vel(1,p) + vel(2,p)*vel(2,p) + &
               vel(3,p)*vel(3,p))
	  mod_a=abs(-2.*ap2*(vel(1,p)*vel(1,p) + vel(2,p)*vel(2,p) + &
               vel(3,p)*vel(3,p)) + bp2*(acc(1,p)*acc(1,p) + &
               acc(2,p)*acc(2,p)+acc(3,p)*acc(3,p)))
	  IF(mod_a.ge.max_acc_local)THEN
	   max_acc_local=mod_a
	   ibod=p
	  ENDIF
	  IF(mod_acc.ge.max_mod_acc)THEN 
	  max_mod_acc=mod_acc
	  ibod2=p
	  ENDIF
	  IF(mod_vel.ge.max_mod_vel) THEN
	  max_mod_vel=mod_vel
	  ibod3=p
	  ENDIF
        ENDDO
    
	max_acc_local=sqrt(max_acc_local)
	
	
	NLONG=1
       CALL MPI_ALLREDUCE(max_acc_local, max_dt_acc, NLONG, MPI_REAL8,  &
        MPI_MAX, MPI_COMM_WORLD, ierror)		
	
       CALL MPI_ALLREDUCE(max_mod_acc, max_modg_acc, NLONG, MPI_REAL8,  &
        MPI_MAX, MPI_COMM_WORLD, ierror)	
       CALL MPI_ALLREDUCE(max_mod_vel, max_modg_vel, NLONG, MPI_REAL8,  &
        MPI_MAX, MPI_COMM_WORLD, ierror)	
		
	   dt_guess=sqrt(0.5*eps/max_dt_acc)
	
	ELSE  ! if dt_var
	
	 dt_guess=dtime	

	ENDIF ! if dt_var
!-----------------------------------------------------------------------------
!   znow determination for output
!-----------------------------------------------------------------------------
	wr_out32=.FALSE.

! We now use the new dt
	
	t_new=tnow+dt_guess
	a_guess=(t_new)**(1.0/alpha)
	z_new=1.0/a_guess-1.0
	
	dt_act=dt_guess
        
	IF(z_new.le.z32_ou(next_out)) THEN   
	
		dt_act=(1/(1+z32_ou(next_out)))**alpha - tnow
		IF(dt_act.eq.0.0) dt_act=dt_guess
!
		wr_out32=.TRUE.
				
		IF(next_out.eq.z32_last)THEN
		   halt_sim=1  !stop simulation at the end  of this time-step
	           IF(me.eq.0) THEN
	 write(uterm,*)'SET STOP SIMULATION: at the end of the current time-step'
	             call flush(uterm)
	           ENDIF 
		ENDIF  
	ENDIF
	
	dtime=dt_act
	dtime2=dtime/2
	
        RETURN
        END