1   subroutine foo(func,p,eval)
2      real(kind=kind(1.0d0)), dimension(3,0:4,0:4,0:4) :: p
3      logical(kind=kind(.true.)), dimension(5,5,5) :: eval
4      interface
5         subroutine func(values,pt)
6            real(kind=kind(1.0d0)), dimension(:), intent(out) :: values
7            real(kind=kind(1.0d0)), dimension(:,:), intent(in) :: pt
8         end subroutine
9      end interface
10      real(kind=kind(1.0d0)), dimension(125,3) :: pt
11      integer(kind=kind(1)) :: n_pt
12
13      n_pt = 1
14      pt(1:n_pt,:) = &
15         reshape( &
16            pack( &
17               transpose(reshape(p,(/3,125/))), &
18               spread(reshape(eval,(/125/)),dim=2,ncopies=3)), &
19            (/n_pt,3/))
20
21   end subroutine
22   end
23