reblock_data Subroutine

public subroutine reblock_data(this, blocklength)

Arguments

Type IntentOptional Attributes Name
real(kind=dp), allocatable :: this(:)
integer, intent(in) :: blocklength

Contents

Source Code


Source Code

    subroutine reblock_data(this, blocklength)
        ! Reduces the data vector for Flyvbjerg and Petersen blocking analysis
        ! by a factor of blocklength (commonly 2, only tested for 2)
        ! General routine, does not require global data

        real(dp), allocatable :: this(:)
        integer :: length, new_length, ind_end, i, j
        real(dp), allocatable :: tmp(:)
        integer, intent(in) :: blocklength
        integer :: ierr
        character(len=*), parameter :: t_r = 'reblock_data'

        if (.not. allocated(this)) call stop_all(t_r, "Error, array not allocated on entry into "&
            & //"reblock_data")
        length = size(this, 1)
        new_length = length / blocklength ! truncates towards zero deliberate, loses data
        !write(stdout,*) "Length is", length
        !write(stdout,*) "Will be reduced to", new_length
        !write(stdout,*) "Loss of data", length-new_length*blocklength
        allocate(tmp(new_length), stat=ierr)
        if (ierr < 0) &
            call stop_all(t_r, 'Bad allocation')
        tmp = 0.0_dp
        j = 1 ! lazy but foolproof - counting elements
        do i = 1, length, blocklength
            ind_end = i + blocklength - 1 ! integer addition is disgusting
            if (ind_end <= length) then
                tmp(j) = average_vector(this(i:ind_end))
            end if
            j = j + 1
        end do
        if (abs(tmp(new_length)) < 1.0e-10_dp) then
            write(stdout, *) 'WARNING: ', t_r, "Whole length of new vector not properly used"
        end if
        deallocate(this)
        allocate(this(new_length))
        this = 0.0_dp
        this = tmp
        deallocate(tmp)

    end subroutine reblock_data