!****************************************************************************
!
!  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, i, j
    real ( kind = 8 ):: TempTimeAveParEne
    real ( kind = 8 ), dimension(:), allocatable :: TempTimeAveEneValley, TempTimeFreqValley, TempTimeRelFreqValley, TempTimeAveParVel
    real ( kind = 8 ), dimension(:), allocatable :: TempTimeRelEneLevelHist
    real ( kind = 8 ), dimension(:,:), allocatable :: TempTimeAveParValleyVel
    real ( kind = 8 ), dimension(:,:), allocatable :: TempTimeRelEneLevelValleyHist
    
    time = 0. ! in simulation module
    step = 0
    flag_write = 0 
    call initialize ()
    call CalculateProperties()
    call WriteAllProperties(step)
    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) )
    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
    do while(time.le.time_max)
          step = step + 1          
          time = time_step*real(step,kind = 8)                 
          call free_flight_scatter()
          ! properties calculation and fileoutput 
          if (time.le.time_early) then  ! for early time
            if (mod(step,step_output1) == 0) then
                flag_write = 1
            else
                flag_write = 0
            endif
          else   ! later time 
            if (mod(step,step_output2) == 0) then
                flag_write = 1
            else
                flag_write = 0
            endif
          endif
          call CalculateProperties() ! TotParEne, AveParEne, TimeTotParEne, TimeAveParEne ! call EveryStepProperties()
          
          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)
            call WriteAllProperties(step) 
            count_step = 0
            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
          else 
            TempTimeAveParEne = TempTimeAveParEne + AveParEne
            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
          endif            
    enddo	! End of the time loop      
end program MonteCarlo