!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!       PERFORM THE FREE-FLIGHT AND SCATTER PART WITHIN ONE TIME INTERVAL   
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine free_flight_scatter()
! common &/ran_var/iso &/scatt_par/emax,de,w(10,3),tau_max(3),max_scatt_mech(3) &/time_1/dt,dtau,tot_time
! &/variables/p(20000,7),ip(20000),energy(20000) &/particle_atr/kx,ky,kz,iv,e &/frequency/freq(10,3)
    use particles
    use scattering
    use simulation
    integer valley,i,j
    real ( kind = 8 ), dimension(DIM) :: Temp_wvk ! carrier position and wave vector k 
    integer :: Temp_valley_i !, Temp_valley_f
    real ( kind = 8 ) :: Temp_ene, Temp_tau
    real ( kind = 8 ):: dte1,dte2,dt2,dt3, dt_remain, err1
 !     Reset counter for scattering frequency
    do i = 1,num_ScatMech 
        do j = 1,num_Valley 
         ScatMechFreq(i,j) = 0.
        enddo
    enddo
    do i = 1, num_Carriers 	! loop for all carriers
!     Inverse mapping of particle atributes
        do j = 1,DIM !col_len = RCutoff_SiIn
            Temp_wvk(j) = Par_wvk(i,j)
        enddo
        Temp_tau = Par_tau(i)
        Temp_valley_i = Par_valley(i)
        Temp_ene = Par_ene(i)        
        if (Temp_ene.ge.1) then
            err1 = 1.0
        endif
!     Initial free-flight of the carriers
        dte1 = Temp_tau
        if(dte1.ge.time_step)then
         dt2= time_step 
        else
         dt2=dte1
        endif      
        call drift(dt2,Temp_wvk,Temp_ene,Temp_valley_i)       
        if(dte1.gt.time_step)goto 401
!     Free-flight and scatter part
402     dte2=dte1        
        call scatter_carrier(Temp_wvk,Temp_ene,Temp_valley_i)          	
219     call random_number(harvest = RandomNum)
        if(RandomNum.le.1e-6) go to 219
        dt3=-(log(RandomNum))*tau_max(Temp_valley_i)
        dt_remain = time_step - dte2	! remaining time to scatter in dt-interval
        if(dt3.le.dt_remain)then
        dt2 = dt3
        else
        dt2 = dt_remain
        endif 
        call drift(dt2,Temp_wvk,Temp_ene,Temp_valley_i)
!     Update times
        dte2 = dte2 + dt3
        dte1 = dte2
        if(dte1.lt.time_step) goto 402
        
401     dte1 = dte1 - time_step 
        Temp_tau = dte1
!     Map particle atributes
        do j = 1,DIM !col_len = RCutoff_SiIn
            Par_wvk(i,j) = Temp_wvk(j)
        enddo
        Par_tau(i) = Temp_tau
        Par_valley(i) = Temp_valley_i
        Par_ene(i) = Temp_ene	             
    enddo
    return	
end