!****************************************************************************
!  PROGRAM: Monte Carlo
!!  PURPOSE:  Entry point for the console application.
!!**************************************************************************** 
program MonteCarlo
    use particles
    use scattering
    use simulation
    implicit none
    integer step,flag_write, count_step, count_bcCheck, count_Poisson, i, j, err1, TempUpdate_Check ! add Tempcheck 11/2012
    real ( kind = 8 ):: TempTimeAveParEne, Time_step_output 
    real ( kind = 8 ), dimension(:), allocatable :: TempTimeAveEneValley, TempTimeFreqValley, TempTimeRelFreqValley, TempTimeAveParVel, TempTimeRelEneLevelHist
    real ( kind = 8 ), dimension(:), allocatable :: Temp_bin_ParNumDen, Temp_bin_ParAveEne, Temp_bin_ParDrftVel
    real ( kind = 8 ), dimension(:), allocatable :: Temp_Node_Potential, Temp_Node_PoiPot, Temp_Node_Charge, Temp_Node_EField ! Poisson
    real(kind=8),dimension(:),allocatable::Temp_Bin_TemperAP,Temp_Bin_TemperOP,Temp_Bin_QsupA,Temp_Bin_QsupO,Temp_Bin_QsupTotal,Temp_Bin_epWm2,Temp_Bin_ppWm2 ! phonon simulation 11/2012 
    real ( kind = 8 ), dimension(:,:), allocatable :: TempTimeAveParValleyVel, TempTimeRelEneLevelValleyHist
    real ( kind = 8 ), dimension(:,:), allocatable :: Temp_binhist_ParNumDen, Temp_bin_ParVal
    integer, dimension(:), allocatable :: Temp_bin_ParNumber
    integer, dimension(:,:), allocatable :: Temp_binhist_ParNumber
    character*30 out_file_1
    call initialize ()
!    call initial_PhononTdis()
    time = 0. ! in simulation module
    step = 0 
    flag_write = 0     
    allocate ( TempTimeAveEneValley(num_Valley) )
    allocate ( TempTimeFreqValley(num_Valley) )
    allocate ( TempTimeRelFreqValley(num_Valley) )
    allocate ( TempTimeAveParVel(num_Valley) )  
    allocate ( TempTimeRelEneLevelHist(num_EnergyLevel) )  
    allocate ( TempTimeAveParValleyVel(num_Valley,DIM) )  
    allocate ( TempTimeRelEneLevelValleyHist(num_Valley,num_EnergyLevel) )
    allocate ( Temp_bin_ParNumDen(num_Bins) )
    allocate ( Temp_bin_ParAveEne(num_Bins) )
    allocate ( Temp_bin_ParVal(num_Valley,num_Bins) )
    allocate ( Temp_bin_ParNumber(num_Bins) ) 
    allocate ( Temp_bin_ParDrftVel(num_Bins) ) 
    allocate ( Temp_Node_Potential(num_Nodes) ) ! Poisson
    allocate ( Temp_Node_PoiPot(num_Nodes) ) 
    allocate ( Temp_Node_Charge(num_Nodes) ) 
    allocate ( Temp_Node_EField(num_Nodes) ) 
    allocate ( Temp_binhist_ParNumDen(num_Bins,num_EnergyLevel) )
    allocate ( Temp_binhist_ParNumber(num_Bins,num_EnergyLevel) )
    count_step = 0
    count_bcCheck = 0
    count_Poisson = 0 
    TempUpdate_Check = 0 ! add Tempcheck 11/2012        
    TempTimeAveParEne = 0.0d0
    TempTimeAveEneValley(1:num_Valley) = 0.0d0
    TempTimeFreqValley(1:num_Valley) = 0.0d0
    TempTimeRelFreqValley(1:num_Valley) = 0.0d0
    TempTimeAveParVel(1:DIM) = 0.0d0
    TempTimeAveParValleyVel(1:num_Valley,1:DIM) = 0.0d0
    TempTimeRelEneLevelHist(1:num_EnergyLevel) = 0.0d0
    TempTimeRelEneLevelValleyHist(1:num_Valley,1:num_EnergyLevel) = 0.0d0
    Temp_bin_ParNumDen(1:num_Bins) = 0.0d0
    Temp_bin_ParAveEne(1:num_Bins) = 0.0d0
    Temp_bin_ParDrftVel(1:num_Bins) = 0.0d0
    Temp_bin_ParVal(1:num_Valley,1:num_Bins) = 0.0d0
    Temp_bin_ParNumber(1:num_Bins) = 0
    Temp_Node_Potential(1:num_Bins) = 0.0d0 !Poisson  
    Temp_Node_PoiPot(1:num_Bins) = 0.0d0
    Temp_Node_Charge(1:num_Bins) = 0.0d0
    Temp_Node_EField(1:num_Bins) = 0.0d0
    Temp_binhist_ParNumDen(1:num_Bins,1:num_EnergyLevel) = 0.0d0
    Temp_binhist_ParNumber(1:num_Bins,1:num_EnergyLevel) = 0     !  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 err1 = 1  ! for debugging endif  endif  enddo 
    ! phonon simulation 11/2012
    allocate ( Temp_Bin_TemperAP(num_Bins) ) 
    allocate ( Temp_Bin_TemperOP(num_Bins) ) 
    allocate ( Temp_Bin_QsupA(num_Bins) ) 
    allocate ( Temp_Bin_QsupO(num_Bins) ) 
    allocate ( Temp_Bin_QsupTotal(num_Bins) ) 
    allocate ( Temp_Bin_epWm2(num_Bins) ) 
    allocate ( Temp_Bin_ppWm2(num_Bins) ) 
    Temp_Bin_TemperAP(1:num_Bins) = 0.0d0 
    Temp_Bin_TemperOP(1:num_Bins) = 0.0d0 
    Temp_Bin_QsupA(1:num_Bins) = 0.0d0
    Temp_Bin_QsupO(1:num_Bins) = 0.0d0
    Temp_Bin_QsupTotal(1:num_Bins) = 0.0d0
    Temp_Bin_epWm2(1:num_Bins) = 0.0d0
    Temp_Bin_ppWm2(1:num_Bins) = 0.0d0
     
    call InitGeneralParas()
    call CalculateProperties()
    call WriteAllProperties(step)    
    ! call PreLoopForField(TimePreFMax,eeHPABPre2,ave_speed2,ave_jeA2,out_file_1)   ! no current loss condition 
    ! call NoCurrentLossCond(eeHPABPre1,eeHPABPre2,ave_speed1,ave_speed2,ave_jeA1,ave_jeA2,)  ! find No current loss using interpolation
    if (PhononSim.eq.3) then ! Middle temperature  1 2 .. external field adjust       
      call PreLoopSdot(TimePreLMax)     ! Find Sdot 
      call InitGeneralParas()
      call FindQuasiSteadyTdist()     ! find quasi steady state Q and T p,A 
      call InitGeneralParas()
    endif 
  
    do while(time.le.time_max)
      step = step + 1          
      time = time_step*real(step,kind = 8)    
      StepParOutL = 0 ! check Particles out
      StepParOutR = 0   !  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 err1 = 1  ! for debugging endif  endif  enddo 
      if (PhononSim.eq.0) then ! Adding Phonon simulation 11/2012
        call free_flight_scatter()    ! Normal
      else
        TempUpdate_Check = TempUpdate_Check + 1 ! Adding Phonon simulation 11/2012
        call free_flight_scatterVar() ! VarPhonon 
        if (TempUpdate_Check.ge.TempUpdateFreq) then 
          TempUpdate_Check = 0 ! add Tempcheck 11/2012
          call Update_Tdist()
          TempVar_BinScatCount(:,:,:) = 0  ! reset 
        endif  
      endif 
     ! properties calculation and fileoutput 
      if (time.le.time_early) then  ! for early time
        Time_step_output = time_step*real(step_output1,kind = 8) 
        if (mod(step,step_output1) == 0) then
          flag_write = 1
        else
          flag_write = 0
        endif
      else   ! later time 
        Time_step_output = time_step*real(step_output2,kind = 8) 
        if (mod(step,step_output2) == 0) then
          flag_write = 1
        else
          flag_write = 0
        endif
      endif
      ! for boundary treat... 
      count_bcCheck = count_bcCheck + 1
      if (Algo_bcCheck.eq.2) then 
          AccumInjParL = AccumInjParL + InjParPerStep + StepParOutL ! elseif (Algo_bcCheck.eq.3) then  AccumInjPar = AccumInjPar + InjParPerStep
          AccumInjParR = AccumInjParR - InjParPerStep + StepParOutR
      elseif (Algo_bcCheck.eq.3) then 
          AccumInjParL = AccumInjParL + num_InjLB ! elseif (Algo_bcCheck.eq.3) then  AccumInjPar = AccumInjPar + InjParPerStep
          AccumInjParR = AccumInjParR + num_InjRB
      endif 
      if (Algo_bcCheck.ne.4) then ! periodic no boundary
        if ((BCoption.eq.2).or.(BCoption.eq.3)) then
          if (count_bcCheck.ge.num_CheckFreq) then
            count_bcCheck = 0
            call BoundaryTreat() 
          endif
        elseif ((BCoption.eq.6).or.(BCoption.eq.7)) then  
          if (count_bcCheck.ge.num_CheckFreq) then
            count_bcCheck = 0
            call BoundaryTreat() 
          endif
        endif 
      else 
        count_bcCheck = 0
      endif 
      ! Poisson Field Update -PoissonSetting, Algo_AssignCharge, Step_FreqFieldUpdate, time_Poisson, count_Poisson
      if ((PoissonSetting.eq.1).and.(time.ge.time_Poisson)) then
        count_Poisson = count_Poisson + 1
        if (count_Poisson.ge.Step_FreqFieldUpdate) then
          count_Poisson = 0
          call ChargeAssign()  ! 
          call PotFromPoisson() !
        endif 
      endif          
      if (time.le.time_early) then 
          AccumParInNumBoundL = 0
          AccumParInNumBoundR = 0
          AccumParOutNumBoundL = 0
          AccumParOutNumBoundR = 0
          AccumParNumBoundPlus = 0
          AccumEneBoundL = 0.0d0  
          AccumEneBoundR = 0.0d0  
          AccumEneBound  = 0.0d0 
      endif !  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 err1 = 1  ! for debugging endif  endif  enddo 
      
      call CalculateProperties() ! !  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 err1 = 1  ! for debugging endif  endif  enddo 
      if (Algo_bcCheck.ne.4) then ! periodic no boundary
        if ((BCoption.eq.0).or.(BCoption.eq.1)) then
          if (count_bcCheck.ge.num_CheckFreq) then
            count_bcCheck = 0
          call BoundaryTreat() 
          endif
        elseif ((BCoption.eq.4).or.(BCoption.eq.5)) then  
          if (count_bcCheck.ge.num_CheckFreq) then
            count_bcCheck = 0
            call BoundaryTreat() 
          endif
        endif
      else 
        count_bcCheck = 0
      endif       
      ! Poisson.. Node properties
      if (PoissonSetting.eq.1) then
        do i = 1,num_Nodes
          Temp_Node_Potential(i) = Temp_Node_Potential(i) + InitPot_Node(i) + pot_Node(i) 
          Temp_Node_PoiPot(i) = Temp_Node_PoiPot(i) + pot_Node(i) 
          Temp_Node_Charge(i) = Temp_Node_Charge(i) + charge_Node(i) 
          Temp_Node_EField(i) = Temp_Node_EField(i) + field_Node(i) 
        enddo
      endif
      ! Phonon simulation 11/2012
      if (PhononSim.ne.0) then
        do i = 1,num_Bins
          Temp_Bin_TemperAP(i) = Temp_Bin_TemperAP(i) + Bin_TemperAP(i)
          Temp_Bin_TemperOP(i) = Temp_Bin_TemperOP(i) + Bin_TemperOP(i)
          Temp_Bin_QsupA(i) = Temp_Bin_QsupA(i) + BinQsupA(i)
          Temp_Bin_QsupO(i) = Temp_Bin_QsupO(i) + BinQsupO(i)
          Temp_Bin_QsupTotal(i) = Temp_Bin_QsupTotal(i) + BinQsupTotal(i) 
          Temp_Bin_epWm2(i) = Temp_Bin_epWm2(i) + epInterEneWm2(i) 
          Temp_Bin_ppWm2(i) = Temp_Bin_ppWm2(i) + ppInterEneWm2(i) 
        enddo 
      endif
      
      TempTimeAveParEne = TempTimeAveParEne + AveParEne
      do i = 1,num_Bins
        Temp_bin_ParNumber(i) = Temp_bin_ParNumber(i) + bin_ParNumber(i)
        Temp_bin_ParNumDen(i) = Temp_bin_ParNumDen(i) + bin_ParNumDen(i)
        Temp_bin_ParAveEne(i) = Temp_bin_ParAveEne(i) + bin_ParAveEne(i)
        Temp_bin_ParDrftVel(i) = Temp_bin_ParDrftVel(i) + bin_ParDrftVel(i)
        do j = 1, num_Valley
          Temp_bin_ParVal(j,i) = Temp_bin_ParVal(j,i) + bin_ParVal(j,i)
        enddo
        do j = 1,num_EnergyLevel
          Temp_binhist_ParNumDen(i,j) = Temp_binhist_ParNumDen(i,j) + binhist_ParNumDen(i,j)
          Temp_binhist_ParNumber(i,j) = Temp_binhist_ParNumber(i,j) + binhist_ParNumber(i,j)
        enddo
      enddo
      do i = 1, num_Valley
        TempTimeAveEneValley(i) = TempTimeAveEneValley(i) + AveEneValley(i)
        TempTimeFreqValley(i) = TempTimeFreqValley(i) + real(FreqValley(i),kind = 8)
        TempTimeRelFreqValley(i) = TempTimeRelFreqValley(i) + RelFreqValley(i)
        TempTimeAveParVel(i) = TempTimeAveParVel(i) + AveParVel(i)
        do j = 1,DIM
          TempTimeAveParValleyVel(i,j) = TempTimeAveParValleyVel(i,j) + AveParValleyVel(i,j)
        enddo
        do j = 1,num_EnergyLevel
          TempTimeRelEneLevelValleyHist(i,j) = TempTimeRelEneLevelValleyHist(i,j) + RelEneLevelValleyHist(i,j)
        enddo
      enddo
      do j = 1,num_EnergyLevel
        TempTimeRelEneLevelHist(j) = TempTimeRelEneLevelHist(j) + RelEneLevelHist(j)
      enddo         
      count_step = count_step + 1
      if (flag_write == 1) then 
          ! TimeAverage Caculation 
          TimeAveParEne = TempTimeAveParEne/real(count_step,kind = 8)
          TimeAveEneValley(1:num_Valley) = TempTimeAveEneValley(1:num_Valley)/real(count_step,kind = 8)
          TimeFreqValley(1:num_Valley) = TempTimeFreqValley(1:num_Valley)/real(count_step,kind = 8)
          TimeRelFreqValley(1:num_Valley) = TempTimeRelFreqValley(1:num_Valley)/real(count_step,kind = 8)
          TimeAveParVel(1:DIM) = TempTimeAveParVel(1:DIM)/real(count_step,kind = 8)
          TimeAveParValleyVel(1:num_Valley,1:DIM) = TempTimeAveParValleyVel(1:num_Valley,1:DIM)/real(count_step,kind = 8)
          TimeRelEneLevelHist(1:num_EnergyLevel) = TempTimeRelEneLevelHist(1:num_EnergyLevel)/real(count_step,kind = 8)
          TimeRelEneLevelValleyHist(1:num_Valley,1:num_EnergyLevel) = TempTimeRelEneLevelValleyHist(1:num_Valley,1:num_EnergyLevel)/real(count_step,kind = 8)            
          TA_bin_ParNumDen(1:num_Bins) = Temp_bin_ParNumDen(1:num_Bins)/real(count_step,kind = 8)
          TA_bin_ParAveEne(1:num_Bins) = Temp_bin_ParAveEne(1:num_Bins)/real(count_step,kind = 8)
          TA_bin_ParVal(1:num_Valley,1:num_Bins) = Temp_bin_ParVal(1:num_Valley,1:num_Bins)/real(count_step,kind = 8)
          TA_bin_ParNumber(1:num_Bins) = real(Temp_bin_ParNumber(1:num_Bins),kind = 8)/real(count_step,kind = 8) 
          TA_bin_ParDrftVel(1:num_Bins) = real(Temp_bin_ParDrftVel(1:num_Bins),kind = 8)/real(count_step,kind = 8) 
          TA_binhist_ParNumDen(1:num_Bins,1:num_EnergyLevel) = Temp_binhist_ParNumDen(1:num_Bins,1:num_EnergyLevel)/real(count_step,kind = 8)
          TA_binhist_ParNumber(1:num_Bins,1:num_EnergyLevel) = real(Temp_binhist_ParNumber(1:num_Bins,1:num_EnergyLevel),kind = 8)/real(count_step,kind = 8)           
          SysCurrent =  real(ParOutNumBoundR,kind = 8)/(count_step*real(step,kind = 8))
          ! Phonon 0709
          if (PrintoutPhonon.eq.1) then ! phonon counting 0709 
            call PhononCountingAndEnergy(Time_step_output)
          endif 
          ! Poisson.. Node properties
          if (PoissonSetting.eq.1) then
            TA_Node_Potential(1:num_Nodes)= Temp_Node_Potential(1:num_Nodes)/real(count_step,kind = 8)
            TA_Node_PoiPot(1:num_Nodes) =  Temp_Node_PoiPot(1:num_Nodes)/real(count_step,kind = 8)
            TA_Node_Charge(1:num_Nodes) =  Temp_Node_Charge(1:num_Nodes)/real(count_step,kind = 8)
            TA_Node_EField(1:num_Nodes) =  Temp_Node_EField(1:num_Nodes)/real(count_step,kind = 8)
          endif
          ! Phonon simulation 11/2012
          if (PhononSim.ne.0) then
            TABin_TemperAP(1:num_Bins) = Temp_Bin_TemperAP(1:num_Bins)/real(count_step,kind = 8)
            TABin_TemperOP(1:num_Bins) = Temp_Bin_TemperOP(1:num_Bins)/real(count_step,kind = 8)
            TABinQsupA(1:num_Bins) = Temp_Bin_QsupA(1:num_Bins)/real(count_step,kind = 8)
            TABinQsupO(1:num_Bins) = Temp_Bin_QsupO(1:num_Bins)/real(count_step,kind = 8)
            TABinQsupTotal(1:num_Bins) = Temp_Bin_QsupTotal(1:num_Bins)/real(count_step,kind = 8)
            TAepInterEneWm2(1:num_Bins) = Temp_Bin_epWm2(1:num_Bins)/real(count_step,kind = 8)
            TAppInterEneWm2(1:num_Bins) = Temp_Bin_ppWm2(1:num_Bins)/real(count_step,kind = 8)
          endif
          if (time.gt.time_early) then
             AccumSysCurrent = real(AccumParOutNumBoundR,kind = 8)/(time-time_early)
          else 
             AccumSysCurrent = 0.0d0
          endif           
          if (time.gt.time_fave) then ! 0720 
             call PropertiesForSummary()
          endif           
          call WriteAllProperties(step) 
          count_step = 0
          ParInNumBoundR = 0
          ParInNumBoundL = 0
          ParOutNumBoundL = 0
          ParOutNumBoundR = 0
          ParNumBoundPlus = 0
          EneBoundL = 0.0d0
          EneBoundR = 0.0d0
          EneBound = 0.0d0
          TempTimeAveParEne = 0.0d0
          TempTimeAveEneValley(1:num_Valley) = 0.0d0
          TempTimeFreqValley(1:num_Valley) = 0.0d0
          TempTimeRelFreqValley(1:num_Valley) = 0.0d0
          TempTimeAveParVel(1:DIM) = 0.0d0
          TempTimeAveParValleyVel(1:num_Valley,1:DIM) = 0.0d0
          TempTimeRelEneLevelHist(1:num_EnergyLevel) = 0.0d0
          TempTimeRelEneLevelValleyHist(1:num_Valley,1:num_EnergyLevel) = 0.0d0
          Temp_bin_ParNumDen(1:num_Bins) = 0.0d0
          Temp_bin_ParAveEne(1:num_Bins) = 0.0d0
          Temp_bin_ParDrftVel(1:num_Bins) = 0.0d0
          Temp_bin_ParVal(1:num_Valley,1:num_Bins) = 0.0d0
          Temp_bin_ParNumber(1:num_Bins) = 0
          Temp_binhist_ParNumDen(1:num_Bins,1:num_EnergyLevel) = 0.0d0
          Temp_binhist_ParNumber(1:num_Bins,1:num_EnergyLevel) = 0
          ! Phonon 0709
          if (PrintoutPhonon.eq.1) then ! phonon counting 0709 
            AllScatCount(:,:) = 0
            BinScatCount(:,:,:) = 0 
          endif          
          ! Poisson.. Node properties
          if (PoissonSetting.eq.1) then
            Temp_Node_Potential(1:num_Nodes) = 0.0d0
            Temp_Node_PoiPot(1:num_Nodes) = 0.0d0
            Temp_Node_Charge(1:num_Nodes) = 0.0d0
            Temp_Node_EField(1:num_Nodes) = 0.0d0
          endif
          ! Phonon simulation 11/2012
          if (PhononSim.ne.0) then
            Temp_Bin_TemperAP(1:num_Bins) = 0.0d0
            Temp_Bin_TemperOP(1:num_Bins) = 0.0d0
            Temp_Bin_QsupA(1:num_Bins) = 0.0d0
            Temp_Bin_QsupO(1:num_Bins) = 0.0d0
            Temp_Bin_QsupTotal(1:num_Bins) = 0.0d0
            Temp_Bin_epWm2(1:num_Bins) = 0.0d0
            Temp_Bin_ppWm2(1:num_Bins) = 0.0d0          
          endif
      endif        
    enddo	! End of the time loop      
end program MonteCarlo