!       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, Temp_pos  ! carrier position and wave vector k 
    integer :: Temp_nAlBin,Temp_valley,Temp_use, Temp_ScatMech, CountValley, bin !, Temp_valley_f 0709
    real ( kind = 8 ) :: Temp_ene, Temp_tau
    real ( kind = 8 ):: timeFreeFlight1,timeFreeFlight2,timeDrift,timeNewFlight,dt_remain,err1
    do i = 1,num_ScatMech  !     Reset counter for scattering frequency
      do j = 1,num_Valley 
        ScatMechFreq(i,j) = 0.
      enddo
    enddo   !    do i = 1, num_maxCarriers if (Par_use(i).eq.1) then   if ((Par_use(i).le.MaxE_eV).and.(Par_use(i).gt.0)) then  else  err1 = 1  ! for debugging endif  endif    enddo
    do i = 1, num_maxCarriers ! num_Carriers 	! loop for all carriers
     if (Par_use(i).eq.1) then !     Inverse mapping of particle atributes
        do j = 1,DIM !col_len = RCutoff_SiIn
            Temp_wvk(j) = Par_wvk(i,j)
            Temp_pos(j) = Par_pos(i,j)
        enddo
        Temp_nAlBin = Par_Al(i)
        Temp_tau = Par_tau(i)
        Temp_valley = Par_valley(i)
        Temp_ene = Par_ene(i)        
        Temp_use = Par_use(i)
        timeFreeFlight1 = Temp_tau !     Initial free-flight of the carriers
        if(timeFreeFlight1.ge.time_step)then
         timeDrift = time_step 
        else
         timeDrift = timeFreeFlight1
        endif 
        call drift(timeDrift,Temp_wvk,Temp_pos,Temp_ene,Temp_valley,Temp_nAlBin,Temp_use) ! Add 08/2013 for Al Content 
        if ((Temp_ene.ge.0).and.(Temp_ene.le.1)) then
            err1 = 0
        else
            err1 = 1
        endif 
        if(Temp_use.eq.0)goto 403       
        if(timeFreeFlight1.gt.time_step)goto 401  !     Free-flight and scatter part
402     timeFreeFlight2=timeFreeFlight1   !        if ((Temp_ene.le.MaxE_eV).and.(Temp_ene.gt.0)) then         else            err1 = 1  ! for debugging        endif
        CountValley = Temp_valley
        call scatter_carrier(Temp_nAlBin,Temp_wvk,Temp_ene,Temp_valley,Temp_ScatMech)    ! 0709 
        if (Temp_ScatMech.ne.0) then ! phonon counting 0709 
            AllScatCount(Temp_ScatMech,CountValley) = AllScatCount(Temp_ScatMech,CountValley) + 1
            bin = num_Bins
            do j = 1,num_Bins ! 
              if ((Temp_pos(1).ge.(loc_bin(j)-0.5*BinSize)).and.(Temp_pos(1).lt.(loc_bin(j)+0.5*BinSize))) then
                bin = j
              endif           
            enddo
            BinScatCount(bin,Temp_ScatMech,CountValley) = BinScatCount(bin,Temp_ScatMech,CountValley) + 1
        endif               	
219     call random_number(harvest = RandomNum)
        if(RandomNum.le.1e-6) goto 219
        timeNewFlight=-(log(RandomNum))*tau_max(Temp_nAlBin,Temp_valley)
        dt_remain = time_step - timeFreeFlight2	! remaining time to scatter in dt-interval
        if(timeNewFlight.le.dt_remain)then
          timeDrift = timeNewFlight
        else
          timeDrift = dt_remain
        endif   !       if ((Temp_ene.le.MaxE_eV).and.(Temp_ene.gt.0)) then    else    err1 = 1  ! for debugging        endif
        call drift(timeDrift,Temp_wvk,Temp_pos,Temp_ene,Temp_valley,Temp_nAlBin,Temp_use) ! Add 08/2013 for Al Content 
        if(Temp_use.eq.0)goto 403 ! if not use... 
!     Update times
        timeFreeFlight2 = timeFreeFlight2 + timeNewFlight
        timeFreeFlight1 = timeFreeFlight2
        if(timeFreeFlight1.lt.time_step) goto 402        
401     timeFreeFlight1 = timeFreeFlight1 - time_step 
        Temp_tau = timeFreeFlight1
!     Map particle atributes
        do j = 1,DIM ! col_len = RCutoff_SiIn
          Par_wvk(i,j) = Temp_wvk(j)
          Par_pos(i,j) = Temp_pos(j)
        enddo
        Par_tau(i) = Temp_tau
        Par_valley(i) = Temp_valley
        Par_Al(i) = Temp_nAlBin
        Par_ene(i) = Temp_ene   ! if ((Par_ene(i).le.MaxE_eV).and.(Par_ene(i).gt.0)) then         else            err1 = 1  ! for debugging        endif
403     Par_use(i) = Temp_use	             
     endif
    enddo
    return	
end
! Al Content Variation Add 08/2013 
subroutine DecideAlBin(posX,nAlBin,LowAlBound,UpAlBound)
  use particles
  use scattering
  use simulation
  real ( kind = 8 ) :: posX,LowAlBound,UpAlBound
  real ( kind = 8 ) :: tempBinMaxX, tempBinMinX
  integer ::  i, nAlBin 
  foundAlBin = 0 
  if (AlContVary.eq.1) then 
    do i = 1,num_BinsForAl  
        tempBinMinX = loc_Start_Al + (i-1)*BinLengthForAl
        if (i.lt.num_BinsForAl) then
            tempBinMaxX = loc_Start_Al + i*BinLengthForAl
        else
            tempBinMaxX = box_len(1) + box_coor(2,1) ! box_coor(2,1) is max value
        endif
        if ((posX.ge.tempBinMinX).and.(posX.lt.tempBinMaxX)) then
            LowAlBound = tempBinMinX
            UpAlBound = tempBinMaxX
            nAlBin = i     
            foundAlBin  = 1  
        endif 
    enddo
  endif 
  if (foundAlBin.eq.0) then 
    nAlBin = num_BinsForAl +1  
    tempBinMinX = box_coor(1,1) - box_len(1)
    if ((AlContVary.eq.1).and.(num_BinsForAl.gt.1)) then 
      tempBinMaxX = loc_Start_Al
    else
      tempBinMaxX = box_len(1) + box_coor(2,1)  
    endif
    foundAlBin  = 1 
    LowAlBound = tempBinMinX
    UpAlBound = tempBinMaxX 
  endif 
  return
end