!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     PERFORM THE K-SPACE AND REAL-SPACE MOTION OF THE CARRIERS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine drift(tau,Temp_wvk,Temp_ene,Temp_valley_i)
    use particles
    use scattering
    use simulation

    real ( kind = 8 ), dimension(DIM) :: Temp_wvk, init_Temp_wvk ! carrier position and wave vector k 
    integer :: Temp_valley_i !, Temp_valley_f
    real ( kind = 8 ) :: Temp_ene, init_Temp_ene !, Temp_tau

    real ( kind = 8 ):: tau
    real ( kind = 8 ):: chahbar, sq_wvk_sum, abs_wvk, factor, err1
    integer i  
    real ( kind = 8 ), dimension(DIM) :: d_wvk, d_pos, sq_wvk

    init_Temp_ene = Temp_ene
    init_Temp_wvk(1) = Temp_wvk(1)
    init_Temp_wvk(2) = Temp_wvk(2)
    init_Temp_wvk(3) = Temp_wvk(3)
    
    !!!!!!!!Test
!    sq_wvk(1) = Temp_wvk(1)*Temp_wvk(1)        
!    sq_wvk(2) = Temp_wvk(2)*Temp_wvk(2)
!    sq_wvk(3) = Temp_wvk(3)*Temp_wvk(3)
!    sq_wvk_sum = sq_wvk(1) + sq_wvk(2) + sq_wvk(3)
!        ! need to consider pos    
 !   abs_wvk = sqrt(sq_wvk_sum)
 !   factor = hhem(Temp_valley_i)*sq_wvk_sum
 !   Temp_ene = 2*factor/(1.+sqrt(1.+nonParabolPara4(Temp_valley_i)*factor))
 !   if (Temp_ene.ge.1) then
 !           err1 = 1.0
 !   endif
    !!!!!!!! Test 
        
    chahbart = echbar*ParCharge*tau
!    Temp_pos, Temp_wvk
!    p(3,i,j)=p(3,i-1,j)-eField(X(3,i-1),F,l,A,ec,bF,barrierCondition)*t; 
!    X(3,i)=X(3,i-1)+p(3,i-1,j)/mee*t-eField(X(3,i-1),F,l,A,ec,bF,barrierCondition)*t^2/2/mee;
    sq_wvk_sum = 0.
    d_wvk(1) = -chahbart*EField(1)
    d_wvk(2) = -chahbart*EField(2)
    d_wvk(3) = -chahbart*EField(3)
         
    Temp_wvk(1) = init_Temp_wvk(1) + d_wvk(1)
    Temp_wvk(2) = init_Temp_wvk(2) + d_wvk(2)
    Temp_wvk(3) = init_Temp_wvk(3) + d_wvk(3)
    
    sq_wvk(1) = Temp_wvk(1)*Temp_wvk(1)        
    sq_wvk(2) = Temp_wvk(2)*Temp_wvk(2)
    sq_wvk(3) = Temp_wvk(3)*Temp_wvk(3)
    sq_wvk_sum = sq_wvk(1) + sq_wvk(2) + sq_wvk(3)
        ! need to consider pos
    
    abs_wvk = sqrt(sq_wvk_sum)
    factor = hhem(Temp_valley_i)*sq_wvk_sum
    Temp_ene = 2*factor/(1.+sqrt(1.+nonParabolPara4(Temp_valley_i)*factor))
!    if (Temp_ene.ge.1) then
!            err1 = 1.0
!    endif
    return
end