Can't compile subroutine with array output with f2py

I have a subroutine that I'm writing in Fortran to be compiled with f2py and the compilation is failing. I won't post the full subroutine here but a MWE is:

SUBROUTINE mwe(Vars, nxc, nyc, vCorr)
IMPLICIT NONE
real(kind=8), dimension(:,:,:,:) :: Vars
integer :: nxc, nyc
integer :: dims(4), nv, nt, nx, ny
real(kind=8), intent(out), allocatable :: vCorr(:,:,:,:)

dims = shape(Vars)
nv=dims(1)
nt=dims(2)
nx=dims(3)
ny=dims(4)
allocate(vCorr(nv, nt, 2*nxc+1, 2*nyc+1))

print*,size(vCorr)
print*,size(Vars)
END SUBROUTINE

This fails with

/tmp/tmpy43di1/src.linux-x86_64-2.7/mwe-f2pywrappers.f:30:31:

       call mwe(vars, nxc, nyc, vcorr)
                               1
Error: Actual argument for ‘vcorr’ must be ALLOCATABLE at (1)

Which apparently means that f2py doesn't accept allocatable output arrays. So I tried to circumvent this problem by passing the shape Vars as an array, so vCorr doesn't have to be allocated, which led me to this code

SUBROUTINE mwe(Vars, nxc, nyc, dims, vCorr)
IMPLICIT NONE
real(kind=8), dimension(:,:,:,:) :: Vars
integer :: nxc, nyc
integer :: dims(4)
real(kind=8) :: vCorr(dims(1),dims(2),2*nxc+1,2*nyc+1)

print*,size(vCorr)
print*,size(Vars)
END SUBROUTINE

Which fails with this error

/tmp/tmp0Y1S9x/src.linux-x86_64-2.7/mwemodule.c:296:39: error: called object ‘dims’ is not a function or function pointer
   vcorr_Dims[0]=dims(1),vcorr_Dims[1]=dims(2),vcorr_Dims[2]=2 * nxc + 1,vcorr_Dims[3]=2 * nyc + 1;

After some look around I came across this page which leads me to believe (even though I'm using f2py2, and not 3) that this is a bug.

Any suggestions around this?

Solution

In the first example f2py does not support allocatable array arguments. They do not fit well with Python arrays.

In the other example f2py does not understand dims(1),dims(2) in vCorr(dims(1),dims(2). It does not support array elements there. This is a bug.

As a workaround use scalar variables for the dimensions

SUBROUTINE mwe(Vars, nxc, nyc, dim1, dim2, vCorr)
  integer, parameter :: dp = kind(1.d0)
  real(dp), dimension(:,:,:,:) :: Vars
  integer :: nxc, nyc
  integer :: dim1, dim2
  real(dp) :: vCorr(dim1,dim2,2*nxc+1,2*nyc+1)

Note: kind=8 is ugly and not portable. The true meaning is NOT 8 bytes, although it does correspond to 8 byte reals in many compilers. But not in all of them. Even the good old double precision is better.