program ensemble

  implicit none

  integer :: i, j
  integer :: k, l
  integer :: a
  integer :: hoge
  integer, parameter :: kf=21    !max-wave-number
  integer, parameter :: m=100    !max-time
  integer, parameter :: u=16     !umsemble-number

  real(8) :: t
  real(8) :: x, y

  real(8) :: E(1:kf)
  real(8) :: en_E(1:kf)

  character(4) :: data_number
  character(1) :: data_a
  character(2) :: data_b
  
  do k=1, m

     do i=1, kf
        E(i)=0.0d0
        en_E(i)=0.0d0
     end do
     
     call counter(k, data_number)

     do j=1, u
        
        if ( j < 10 ) then

           call counter_1(j, data_a)
           
           open(100+j,file='test'//data_a//'_energy_spectrum_'//data_number//'.dat')
           do i=1, kf
              read(100+j,*) hoge, E(i)
              en_E(i)=en_E(i)+E(i)
           end do
           close(100+j)

        else

           call counter_2(j, data_b)

           open(100+j,file='test'//data_b//'_energy_spectrum_'//data_number//'.dat')
           do i=1, kf
              read(100+j,*) hoge, E(i)
              en_E(i)=en_E(i)+E(i)
           end do
           close(100+j)

        end if

     end do
     
     open(200,file='ensemble_energy_spectrum_'//data_number//'.dat')
     do i=1, kf
        write(200,*) i, en_E(i)/dble(u)
     end do
     close(200)
     
 end do

end program ensemble

!----- データ出力 -----
subroutine counter(i, data_number)

  implicit none

  integer, intent(in) :: i
  character(4), intent(out) :: data_number
  integer :: o_4, o_3, o_2, o_1
  character(10) :: i_data
  
  i_data='0123456789'

  o_4=i/1000
  o_3=(i-o_4*1000)/100
  o_2=(i-o_4*1000-o_3*100)/10
  o_1= i-o_4*1000-o_3*100-o_2*10
  data_number=i_data(o_4+1:o_4+1)//&
       i_data(o_3+1:o_3+1)//&
       i_data(o_2+1:o_2+1)//&
       i_data(o_1+1:o_1+1)
  
end subroutine counter

!----- データ出力 -----
subroutine counter_1(i, data_a)

  implicit none

  integer, intent(in) :: i
  character(1), intent(out) :: data_a
  integer :: o_1
  character(10) :: i_data
  
  i_data='0123456789'

  o_1=i
  
  data_a=i_data(o_1+1:o_1+1)
  
end subroutine counter_1

!----- データ出力 -----
subroutine counter_2(i, data_b)

  implicit none

  integer, intent(in) :: i
  character(2), intent(out) :: data_b
  integer :: o_2, o_1
  character(10) :: i_data
  
  i_data='0123456789'

  o_2=i/10
  o_1=i-o_2*10
  data_b=i_data(o_2+1:o_2+1)//&
       i_data(o_1+1:o_1+1)
  
end subroutine counter_2
