program statistics
implicit none
integer :: NPoints = 10000
integer :: NBlocks = 10000
integer :: istat, xsum, i
integer, parameter :: dp = selected_real_kind(15,300)
real(kind=dp) :: Xavg  !Average
real(kind=dp) :: Std   !Standard Deviation
real(kind=dp) :: Sign_Level = 0.95
real(kind=dp) :: xbar=0, Nxvar=0, xvar=0, xstd=0
real(kind=dp), dimension(10000) :: X
real(kind=dp), dimension(10000) :: XBavg

open(unit=100, file="KK.txt", status="old", action="read",iostat=istat)
if (istat .ne. 0) stop "error opening file"

open(unit=200, file = 'err.dat')

read(100,*) X

call blockav( NPoints, NBlocks, X, XBavg, Xavg, Std)


print *, Std

call meanvar(x,xbar,xvar,Nxvar, xstd)


contains

 subroutine meanvar(x,xbar,xvar,Nxvar, xstd)
  integer, parameter :: dp = selected_real_kind(15,300)
  real(kind=dp), intent(out) :: xbar, Nxvar, xvar, xstd
  integer :: i, j
  real(kind=dp) :: xsum
  real(kind=dp), dimension(10000), intent(in) :: x

  do i = 1,10000

      xsum = xsum + X(i)

      xvar = Nxvar/10000

  end do

      xbar = xsum/10000

  do i = 1,10000

      Nxvar = Nxvar + (x(i) - xbar)**2
  end do

      xvar = Nxvar/10000
      xstd = xvar**1/2
end subroutine

subroutine blockav( NPoints, NBlocks, X, XBavg, Xavg, Std)

integer,intent(in) :: NPoints   !Number of values to be averaged
integer,intent(in) :: NBlocks       !Number of blocks to be used
integer, parameter :: dp = selected_real_kind(15,300)
real(kind=dp), dimension(NPoints), intent(in) :: X     !Array for values to be block averaged
real(kind=dp), dimension(NBlocks), intent(out) :: XBavg !Block Averages
real(kind=dp), intent(out) :: Xavg  !Average
real(kind=dp), intent(out) :: Std   !Standard Deviation
real(kind=dp):: Sign_Level = 0.95


!Local variables.
integer::  iblock, ipoint
integer::  FirstPoint, Pointsperblock




Pointsperblock = NPoints / NBlocks	 !Integer division

if (Pointsperblock.eq.0) then
	write (*,*) 'The number of Points per block is too large for the specified number of blocks'
	write (*,*) 'Please specify a smaller number of Points per block, or increase the number of Points'
	stop
endif

FirstPoint = MOD( NPoints, NBlocks )


!Initializing
XBavg = 0.0
Xavg = 0.0
Std = 0.0


!Calculate the block averages.
do iblock = 1, NBlocks
	do ipoint = FirstPoint + (iblock - 1) * Pointsperblock + 1, &
	       		FirstPoint + iblock * Pointsperblock
		XBavg(iblock) = XBavg(iblock) + X(ipoint)

	end do
	XBavg(iblock) = XBavg(iblock) / real(Pointsperblock)
	Std = sqrt((Std + ( XBavg(iblock) - Xavg ) ** 2)/(Nblocks -1))
	write(200,*) Std

end do

!Average over the simulation.
Xavg = SUM( XBavg( 1:NBlocks ) ) / real( NBlocks )





return

end subroutine

end program