Is it a bug: data aligment of matrix returned by gemm

This is my first post. Greetings everyone!

I am not sure if the following behaviour is a bug or simply the matter
of my poor understanding of F90 standard so I am not submitting this a 
a formal bug report.

Consider following code:

function matmul_func(a,b)
  real(kind=8),intent(in) :: a(:,:), b(:,:)
  real(kind=8) :: matmul_func(size(a,1),size(b,2))

  call dgemm('N','N',size(a,1),size(b,2),size(a,2),1.0d0,a,size(a,1),&
            &b,size(b,1),0.0d0,matmul_func,size(a,1))

end function

program test
  integer, paramater :: range=2
  real(kind=8) :: a(rank,rank), b(rank,rank), c(2,rank,rank)

  interface
    function matmul_func....
                ....
    end function matmul_func
  end interface

  a(1,:)=1.0d0
  a(2,:)=2.0d0

  b(1,:)=(/1.0d0,2.0d0/)
  b(2,:)=(/3.0d0,4.0d0/)

  c(1,:,:)=matmul_func(a,b)

end program test

The result is incorrect in a sense the the resulting
2x2 matrix is put into c coulmnwise as: 

c(1,1,1), c(2,1,1), c(1,2,1), c(2,2,1) 

whereas the syntax suggests that the c(1,:,:) block should be filled,
i.e.:

c(1,1,1), c(1,2,1), c(1,1,2), c(1,2,2).

Is it a bug? The result can be understood in the above fashion and
corrected for, but it would be nice for gfortran to
honour the block syntax. Intel Fortran and g95 are doing just that.

If I pass the rank of the matrices explicitely, i.e.

 function matmul_func(n,a,b)
	real... :: a(n,n), b(n,n), matmul_func(n,n)

with the rest of the body of the function as previously the results are
correct.

Replacing  geem with: matmul_func=matmul(a,b) yields correct
results. Also when GEMM  is called in the main body of the program:

call dgemm(..... c(1,:,:),rank)

the result is correct.

I have used 4.3,4.4 and the 4.5-20090618 versions linking agains
reference netlib BLAS (also recompiled with 4.5) and intel's mkl and
ATLAS. The results are all the same.

I'd appreciate your comments
Regards
MZ