PrintBlocking Subroutine

    SUBROUTINE PrintBlocking(Iter)
        use util_mod, only: get_free_unit
        INTEGER :: i, NoBlocks, Iter, NoBlockSizes, NoContrib, iunit
        real(dp) :: MeanEn, MeanEnSqrd, StandardDev, Error, ErrorinError
        CHARACTER(len=30) :: abstr

!First find out how many blocks would have been formed with the number of iterations actually performed.

        NoContrib = 0
        IF (tBlockEveryIteration) THEN
            NoContrib = (Iter - StartBlockIter) + 1
        ELSE
            NoContrib = ((Iter - StartBlockIter) / StepsSft) + 1
        end if
        IF (NoContrib == 1) RETURN

        NoBlockSizes = 0
        NoBlockSizes = FLOOR((LOG10(REAL(NoContrib - 1, dp))) / (LOG10(2.0_dp)))

        iunit = get_free_unit()
        IF (tSaveBlocking) THEN
!We are saving the blocking files - write to new file.
            abstr = 'BLOCKINGANALYSIS-'//str(Iter)
            open(iunit, file=abstr, status='unknown')
        ELSE
            open(iunit, file='BLOCKINGANALYSIS', status='unknown')
        end if

        write(iunit, '(I4,A,I4)') NoBlockSizes, ' blocks were formed with sizes from 1 to ', (2**(NoBlockSizes))
        write(iunit, '(3A16,5A20)') '1.Block No.', '2.Block Size  ', '3.No. of Blocks', '4.Mean E', &
            '5.Mean E^2', '6.SD', '7.Error', '8.ErrorinError'

        do i = 0, NoBlockSizes - 1

            ! First need to find out how many blocks of this particular size contributed to the final sum in BlockSum.
            ! NoContrib is the total number of contributions to the blocking throughout the simulation.
            NoBlocks = 0
            NoBlocks = FLOOR(REAL(NoContrib / (2**i), dp))
            ! This finds the lowest integer multiple of 2**i (the block size).

            MeanEn = BlockSum(i) / REAL(NoBlocks, dp)
            MeanEnSqrd = BlockSqrdSum(i) / REAL(NoBlocks, dp)

            StandardDev = SQRT(MeanEnSqrd - (MeanEn**2))
            IF (abs(StandardDev) < 1.0e-12_dp) THEN
                Error = 0.0_dp
                ErrorinError = 0.0_dp
            ELSE
                Error = StandardDev / SQRT(REAL(NoBlocks - 1, dp))
                ErrorinError = Error / SQRT(2.0_dp * (REAL(NoBlocks - 1, dp)))
            end if

!This is from the blocking paper, and indicates the error in the blocking error, due to the limited number of blocks available.

            write(iunit, '(3I16,5F20.10)') i, (2**i), NoBlocks, MeanEn, MeanEnSqrd, StandardDev, Error, ErrorinError

        end do

        close(iunit)

    END SUBROUTINE PrintBlocking