shuffle_3Didata Subroutine

  subroutine shuffle_3Didata(data,dim,seed)
    !! Shuffle a 3D array along one dimension
    implicit none

    ! Arguments
    integer, dimension(:,:,:), intent(inout) :: data
    !! 3D array to be shuffled
    integer, intent(in) :: dim
    !! Dimension along which to shuffle
    integer, optional, intent(in) :: seed
    !! Random seed

    ! Local variables
    integer :: istart,seed_size
    !! Start index, seed size
    integer :: i,j,n_data
    !! Loop indices, number of data points
    real(real32) :: r
    !! Random number
    integer, dimension(3) :: idx_s,idx_e,jdx_s,jdx_e
    !! Start and end indices
    integer, dimension(3,2) :: t_size
    !! Temporary size
    integer, allocatable, dimension(:) :: iseed
    !! Random seed
    integer, allocatable, dimension(:,:,:) :: tlist
    !! Temporary list


    ! Set or get random seed
    !---------------------------------------------------------------------------
    call random_seed(size=seed_size)
    allocate(iseed(seed_size))
    if(present(seed))then
       iseed = seed
       call random_seed(put=iseed)
    else
       call random_seed(get=iseed)
    end if


    ! Get the size of the data
    !---------------------------------------------------------------------------
    n_data = size(data,dim=dim)
    do i=1,3
       t_size(i,1) = 1
       jdx_s(i) = 1
       jdx_e(i) = size(data,dim=i)
       idx_s(i) = 1
       idx_e(i) = size(data,dim=i)
       if(i.eq.dim)then
          t_size(i,2) = 1
       else
          t_size(i,2) = size(data,dim=i)
       end if
    end do
    allocate(tlist(t_size(1,2),t_size(2,2),t_size(3,2)))


    ! Shuffle the data
    !---------------------------------------------------------------------------
    istart=1
    do i=1,n_data
       call random_number(r)
       j = istart + floor((n_data+1-istart)*r)
       if(i.eq.j) cycle
       idx_s(dim) = i
       idx_e(dim) = i
       jdx_s(dim) = j
       jdx_e(dim) = j
       tlist(&
            t_size(1,1):t_size(1,2),&
            t_size(2,1):t_size(2,2),&
            t_size(3,1):t_size(3,2)) = data(&
            idx_s(1):idx_e(1),&
            idx_s(2):idx_e(2),&
            idx_s(3):idx_e(3))
       data(&
            idx_s(1):idx_e(1),&
            idx_s(2):idx_e(2),&
            idx_s(3):idx_e(3)) = data(&
            jdx_s(1):jdx_e(1),&
            jdx_s(2):jdx_e(2),&
            jdx_s(3):jdx_e(3))
       data(&
            jdx_s(1):jdx_e(1),&
            jdx_s(2):jdx_e(2),&
            jdx_s(3):jdx_e(3)) = tlist(&
            t_size(1,1):t_size(1,2),&
            t_size(2,1):t_size(2,2),&
            t_size(3,1):t_size(3,2))
    end do

  end subroutine shuffle_3Didata