shuffle_5Drdata_1Drlist Subroutine

  subroutine shuffle_5Drdata_1Drlist(data,label,dim,seed,shuffle_list)
    !! Shuffle a 5D array along one dimension
    implicit none

    ! Arguments
    real(real32), dimension(:,:,:,:,:), intent(inout) :: data
    !! 5D array to be shuffled
    real(real32), dimension(:), intent(inout) :: label
    !! 1D array to be shuffled
    integer, intent(in) :: dim
    !! Dimension along which to shuffle
    integer, optional, intent(in) :: seed
    !! Random seed
    integer, optional, dimension(size(data,dim)), intent(out) :: shuffle_list
    !! Index array

    ! Local variables
    integer :: istart,seed_size
    !! Start index, seed size
    integer :: i,j,n_data
    !! Loop indices, number of data points
    real(real32) :: rtmp1
    !! Temporary real
    real(real32) :: r
    !! Random number
    integer, dimension(5) :: idx_s,idx_e,jdx_s,jdx_e
    !! Start and end indices
    integer, dimension(5,2) :: t_size
    !! Temporary size
    integer, allocatable, dimension(:) :: iseed
    !! Random seed
    real(real32), 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,5
       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),t_size(4,2),&
         t_size(5,2)))

    istart=1
    do i=1,n_data
       call random_number(r)
       j = istart + floor((n_data+1-istart)*r)
       if(present(shuffle_list)) shuffle_list(i) = j
       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),&
            t_size(4,1):t_size(4,2),&
            t_size(5,1):t_size(5,2)) = data(&
            idx_s(1):idx_e(1),&
            idx_s(2):idx_e(2),&
            idx_s(3):idx_e(3),&
            idx_s(4):idx_e(4),&
            idx_s(5):idx_e(5))
       data(&
            idx_s(1):idx_e(1),&
            idx_s(2):idx_e(2),&
            idx_s(3):idx_e(3),&
            idx_s(4):idx_e(4),&
            idx_s(5):idx_e(5)) = data(&
            jdx_s(1):jdx_e(1),&
            jdx_s(2):jdx_e(2),&
            jdx_s(3):jdx_e(3),&
            jdx_s(4):jdx_e(4),&
            jdx_s(5):jdx_e(5))
       data(&
            jdx_s(1):jdx_e(1),&
            jdx_s(2):jdx_e(2),&
            jdx_s(3):jdx_e(3),&
            jdx_s(4):jdx_e(4),&
            jdx_s(5):jdx_e(5)) = 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),&
            t_size(4,1):t_size(4,2),&
            t_size(5,1):t_size(5,2))
       rtmp1 = label(i)
       label(i) = label(j)
       label(j) = rtmp1
    end do

  end subroutine shuffle_5Drdata_1Drlist