shuffle_2Drdata Subroutine

  subroutine shuffle_2Drdata(data,dim,seed)
    implicit none

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

    ! Local variables
    integer :: istart,seed_size
    !! Start index, seed size
    integer :: i,j,n_data,iother
    !! Loop indices, number of data points, other dimension
    integer :: i1s,i2s,i1e,i2e,j1s,j2s,j1e,j2e
    !! Start and end indices
    real(real32) :: r
    !! Random number
    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


    ! Shuffle the data
    !---------------------------------------------------------------------------
    n_data = size(data,dim=dim)
    if(dim.eq.1)then
       iother = 2
       i2s=1;i2e=size(data,dim=iother)
       j2s=1;j2e=size(data,dim=iother)
    else
       iother = 1
       i1s=1;i1e=size(data,dim=iother)
       j1s=1;j1e=size(data,dim=iother)
    end if
    istart=1
    allocate(tlist(1,size(data,dim=iother)))
    do i=1,n_data
       call random_number(r)
       j = istart + floor((n_data+1-istart)*r)
       if(i.eq.j) cycle
       if(dim.eq.1)then
          i1s=i;i1e=i
          j1s=j;j1e=j
       else
          i2s=i;i2e=i
          j2s=j;j2e=j
       end if
       tlist(1:1,:) = reshape(data(i1s:i1e,i2s:i2e),shape=shape(tlist))
       data(i1s:i1e,i2s:i2e) = data(j1s:j1e,j2s:j2e)
       data(j1s:j1e,j2s:j2e) = reshape(tlist(1:1,:),&
            shape=shape(data(j1s:j1e,j2s:j2e)))
    end do

  end subroutine shuffle_2Drdata