subroutine initialize ( )
!*****************************************************************************
!! INITIALIZE is the initialization procedure.
!*****************************************************************************
  use simulation
  use scattering
  use particles
  use nodes
  implicit none
!  Read the input sample containing the initial particle coordinates and parameters.
  call read_sample()
  call sc_table()
  call initial_printout()
  call initial_distribution()
!  Print information on the run on the standard output file.
  return
end

subroutine read_sample () ! parameter setting 
  use simulation
  use scattering
  use particles
  implicit none
  real ( kind = 8 ) :: int_temp
  character ( len = 80 ) :: temp
  character ( len = 80 ) :: Jobname
  integer :: i, j, p, q
!  Now that we know the system size, we can dynamically allocate the arrays containing atomic informations
  open ( unit = 1, file = 'input.txt', status = 'old', action = 'read', err = 600 )
  read(1,*,end=800,err=900) temp
  read(1,*,end=600,err=900) num_EnergyLevel , num_Carriers, step_output1, step_output2 !, Ntot, NDOSGroup1, NDOSGroup2,TRescaling, box_coor(1:2,1), box_coor(1:2,2), box_coor(1:2,3) 
  read(1,*,end=600,err=900) acoustic_G,polar_G,intervalley_G_L,intervalley_G_X, acoustic_L,polar_L,intervalley_L_G,intervalley_L_L,intervalley_L_X,acoustic_X,polar_X,intervalley_X_G,intervalley_X_L,intervalley_X_X
  read(1,*,end=600,err=900) Te, Tp, Tsys, EField(1:DIM), time_max, time_step, time_early, MaxE_eV, StepE_eV, rho_GaAs, u_sound, ParCharge
  read(1,*,end=600,err=900) rel_eff_mass(1:num_Valley), nonParabolPara(1:num_Valley), n_d, n_i, primV
  read(1,*,end=600,err=900) eps_high,eps_low, acou_sigma_G,acou_sigma_L,acou_sigma_X
  read(1,*,end=600,err=900) DefPot_G_L,DefPot_G_X,DefPot_L_G,DefPot_L_L,DefPot_L_X,DefPot_X_G,DefPot_X_L,DefPot_X_X
  read(1,*,end=600,err=900) polar_Ep_G, polar_Ep_L, polar_Ep_X, IntVal_Ep_G_L, IntVal_Ep_G_X, IntVal_Ep_L_G, IntVal_Ep_L_L, IntVal_Ep_L_X, IntVal_Ep_X_G, IntVal_Ep_X_L, IntVal_Ep_X_X
  read(1,*,end=600,err=900) split_L_G,split_X_G,eq_valleys_G,eq_valleys_L,eq_valleys_X
  read(1,*,end=800,err=900) temp
!  if ( Ntot <= 0 ) then
!    print*,'Read_Sample: FATAL: Ntot is', Ntot
!    stop
!  end if
    allocate ( scatt_table(num_EnergyLevel,num_ScatMech,num_Valley) ) !scatt_table(i,i_count_machanism,valley)
    allocate ( ren_scatt_table(num_EnergyLevel,num_ScatMech,num_Valley) )
    allocate ( Par_wvk(num_Carriers,DIM) )
    allocate ( Par_pos(num_Carriers,DIM) )
    allocate ( Par_vel(num_Carriers,DIM) )
    allocate ( Par_valley(num_Carriers) ) 
    allocate ( Par_ene(num_Carriers) ) 
    allocate ( Par_tau(num_Carriers) ) 
    allocate ( RelCumFreqDist(num_EnergyLevel) ) 
    allocate ( FreqDist(num_EnergyLevel) ) 
    allocate ( CumFreqDist(num_EnergyLevel) )       
    allocate ( TotEneValley(num_Valley) )
    allocate ( AveEneValley(num_Valley) )
    allocate ( TimeAveEneValley(num_Valley) )
    allocate ( FreqValley(num_Valley) )  
    allocate ( RelFreqValley(num_Valley) )  
    allocate ( TimeFreqValley(num_Valley) )  
    allocate ( TimeRelFreqValley(num_Valley) )  
    allocate ( AveParValleyVel(num_Valley,DIM) )    
    allocate ( TimeAveParValleyVel(num_Valley,DIM) )    
    allocate ( EneLevelHist(num_EnergyLevel) )    
    allocate ( RelEneLevelHist(num_EnergyLevel) )    
    allocate ( TimeRelEneLevelHist(num_EnergyLevel) )     
    allocate ( EneLevelValleyHist(num_Valley,num_EnergyLevel) )       
    allocate ( RelEneLevelValleyHist(num_Valley,num_EnergyLevel) )       
    allocate ( TimeRelEneLevelValleyHist(num_Valley,num_EnergyLevel) ) 
    echbar = e_c/hbar 
    eps_high = eps_high*eps_o
    eps_low = eps_low*eps_o
    do i = 1,num_Valley 	
      eff_mass(i) = m_e*rel_eff_mass(i)  
      nonParabolPara2(i) = 2.*nonParabolPara(i)
      nonParabolPara4(i) = 4.*nonParabolPara(i)
    enddo
    do i = 1,num_Valley 	
      semh(i)=sqrt(2.*eff_mass(i))*sqrt(e_c)/hbar
      hhem(i)=hbar/eff_mass(i)/e_c*hbar/2. ! hbar*hbar/(eff_mass(i)*e_c*2.)
    enddo
    Te_eV = Te*kB/e_c
    Tp_eV = Tp*kB/e_c 
    Tsys_eV = Tsys*kB/e_c
! initialize
    max_scat = 0.0;
    flag_mech(:,:) = 0
    Par_wvk(:,:) = 0.0d0 
    Par_pos(:,:) = 0.0d0
    Par_vel(:,:) = 0.0d0
    Par_valley(:) = 0
    Par_ene(:) = 0.0D+00
    Par_tau(:) = 0.0D+00
    scatt_table(:,:,:) = 0.0D+00
    ren_scatt_table(:,:,:) = 0.0D+00
    flag_mech(:,:) = 0
    E_change(:,:) = 0.0d0
    i_valley(:,:) = 0
    ScatMechFreq(:,:) = 0
    tau_max(:) = 0.0d0
    max_scat = 0.0d0
!    Temp_pos(:) = 0.0d0
!    Temp_wvk(:) = 0.0d0
!    Temp_valley_i = 0
!    Temp_valley_f = 0
!    Temp_ene = 0.0d0
!    Temp_tau = 0.0d0
    ParCharge = 1.0;
    ! Fermilevel calculation 
    E_F = kB*Te*log(n_d/n_i)/e_c-0.674;  
    call RandomNumGen() ! Random..   
!  Compute volume and density
! do i = 1, DIM
!     box_len(i) = box_coor(2,i) - box_coor(1,i)
! end do
!  volume = product(box_len(:)) !  density = real ( NTYPE(4), kind = 8 ) / volume
! Temp_o = ep_o/kB
  close ( unit = 1 )
 return
!  Handling of various kinds of errors
600 continue
   print*,'Read_Sample: FATAL: input.txt is empty?'
   stop
700 continue
   print*,'Read_Sample: FATAL: input.txt not found.'
   stop
800 continue
   print*,'Read_Sample: FATAL: premature end-of-file at atom ',i
   close(unit=1)
    stop   
900 continue
   print*,'Read_Sample: FATAL: read error in input.txt'
   close(unit=1)
   stop
end

subroutine initial_printout ( )
implicit none
integer :: i
open ( unit = 2, file = 'TimePointEneFreq.txt', status = 'replace', action = 'write', err = 610 )
  open ( unit = 3, file = 'TimeAveEneFreq.txt', status = 'replace', action = 'write', err = 610 )
  open ( unit = 4, file = 'TimePointVel.txt', status = 'replace', action = 'write', err = 610 )
  open ( unit = 5, file = 'TimeAveVel.txt', status = 'replace', action = 'write', err = 610 )
  open ( unit = 6, file = 'ParticlesForTime.txt', status = 'replace', action = 'write', err = 620 )
  open ( unit = 7, file = 'PointEnergyHist.txt', status = 'replace', action = 'write', err = 620 )
  open ( unit = 15, file = 'TimeAveEnergyHist.txt', status = 'replace', action = 'write', err = 620 )
 ! open ( unit = 14, file = 'qtot.txt', status = 'replace', action = 'write', err = 610 )
 ! open ( unit = 20, file = 'TorFreq.txt', status = 'replace', action = 'write', err = 610 )
 ! open ( unit = 21, file = 'TorTrend.txt', status = 'replace', action = 'write', err = 610 )
 ! open ( unit = 40, file = 'TorFreq2.txt', status = 'replace', action = 'write', err = 610 )
 ! open ( unit = 41, file = 'TorTrend2.txt', status = 'replace', action = 'write', err = 610 )
 ! open ( unit = 50, file = 'TorFreqA.txt', status = 'replace', action = 'write', err = 610 )
 ! open ( unit = 51, file = 'TorTrendA.txt', status = 'replace', action = 'write', err = 610 )
return  
!  Handling of various kinds of errors
610 continue
   print*,'FATAL: Error in Output.xyz.'
   stop
620 continue
   print*,'FATAL: Error in ParticlesForTime.txt'
   stop
630 continue
   print*,'Read_Sample: FATAL: premature end-of-file at atom ',i
   stop   
640 continue
   print*,'Read_Sample: FATAL: read error in input.txt'
   stop
end