Problem with f2py not modifying an inout array

I have a fortran routine that takes 6 arrays as input, and modifies the content of the first array. The arrays are all numpys with the following flags

  C_CONTIGUOUS : False
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

I tried declaring the array that needs to be modified as intent(in, out), also tried intent(inout), and using the inplace modifier, but in every case, once the routine terminates, the array that was supposed to be modified is not modified on the Python side. Using the -DF2PY_REPORT_ON_ARRAY_COPY=1 I also note that 5 arrays are copied. This of course creates a performance issue. What I am doing wrong? This is the fortran subroutine

          subroutine POISSONOP_GSRB(
     & phi
     & ,iphilo0,iphilo1
     & ,iphihi0,iphihi1
     & ,nphicomp
     & ,rhs
     & ,irhslo0,irhslo1
     & ,irhshi0,irhshi1
     & ,nrhscomp
     & ,Mx
     & ,iMxlo0,iMxlo1
     & ,iMxhi0,iMxhi1
     & ,nMxcomp
     & ,My
     & ,iMylo0,iMylo1
     & ,iMyhi0,iMyhi1
     & ,nMycomp
     & ,Mz
     & ,iMzlo0,iMzlo1
     & ,iMzhi0,iMzhi1
     & ,nMzcomp
     & ,Dinv
     & ,iDinvlo0,iDinvlo1
     & ,iDinvhi0,iDinvhi1
     & ,idestBoxlo0,idestBoxlo1
     & ,idestBoxhi0,idestBoxhi1
     & ,whichPass
     & )
      implicit none
      integer nphicomp
      integer iphilo0,iphilo1
      integer iphihi0,iphihi1
      REAL*8 phi(
     & iphilo0:iphihi0,
     & iphilo1:iphihi1,
     & 0:nphicomp-1)     
cf2py intent(in, overwrite) phi
      integer nrhscomp
      integer irhslo0,irhslo1
      integer irhshi0,irhshi1
      REAL*8 rhs(
     & irhslo0:irhshi0,
     & irhslo1:irhshi1,
     & 0:nrhscomp-1)
cf2py intent(in, inplace) rhs
      integer nMxcomp
      integer iMxlo0,iMxlo1
      integer iMxhi0,iMxhi1
      REAL*8 Mx(
     & iMxlo0:iMxhi0,
     & iMxlo1:iMxhi1,
     & 0:nMxcomp-1)
cf2py intent(in, inplace) Mx
      integer nMycomp
      integer iMylo0,iMylo1
      integer iMyhi0,iMyhi1
      REAL*8 My(
     & iMylo0:iMyhi0,
     & iMylo1:iMyhi1,
     & 0:nMycomp-1)
cf2py intent(in, inplace) My
      integer nMzcomp
      integer iMzlo0,iMzlo1
      integer iMzhi0,iMzhi1
      REAL*8 Mz(
     & iMzlo0:iMzhi0,
     & iMzlo1:iMzhi1,
     & 0:nMzcomp-1)
cf2py intent(in, inplace) Mz
      integer iDinvlo0,iDinvlo1
      integer iDinvhi0,iDinvhi1
      REAL*8 Dinv(
     & iDinvlo0:iDinvhi0,
     & iDinvlo1:iDinvhi1)
cf2py intent(in, inplace) Dinv
      integer idestBoxlo0,idestBoxlo1
      integer idestBoxhi0,idestBoxhi1
      integer whichPass
      integer i,j
      integer indtot, imin, imax
      integer n, ncomp
      REAL*8 Sphi
      REAL*8 MyL
      REAL*8 MyR
      ncomp = nphicomp
      do n = 0, ncomp-1
          do j = idestBoxlo1, idestBoxhi1
            MyL = My(0,j,0)
            MyR = My(0,j,1)
            imin = idestBoxlo0
            indtot = imin + j
            imin = imin + abs(mod(indtot + whichPass, 2))
            imax = idestBoxhi0
            do i = imin, imax, 2
              Sphi =
     & Mx(i,0,0) * phi(i-1,j,n)
     & + Mx(i,0,1) * phi(i+1,j,n)
     & + MyL * phi(i,j-1,n)
     & + MyR * phi(i,j+1,n)
              phi(i,j,n) = Dinv(i,j)
     & * (rhs(i,j,n) - Sphi)
      enddo
      enddo
      enddo
      return
      end

and this is how it is called on Python

GS.poissonop_gsrb(phi.data,phi.box.LoEnd()[0],phi.box.LoEnd()[1],phi.box.HiEnd()[0], phi.box.HiEnd()[1],
                      rhs.data,rhs.box.LoEnd()[0],rhs.box.LoEnd()[1],rhs.box.HiEnd()[0], rhs.box.HiEnd()[1],
                      Mx.data,Mx.box.LoEnd()[0],Mx.box.LoEnd()[1],Mx.box.HiEnd()[0], Mx.box.HiEnd()[1],
                      My.data,My.box.LoEnd()[0],My.box.LoEnd()[1],My.box.HiEnd()[0], My.box.HiEnd()[1],
                      Mz.data,Mz.box.LoEnd()[0],Mz.box.LoEnd()[1],Mz.box.HiEnd()[0], Mz.box.HiEnd()[1],
                      DInv.data,DInv.box.LoEnd()[0],DInv.box.LoEnd()[1],DInv.box.HiEnd()[0], DInv.box.HiEnd()[1],
                      L[0], L[1],H[0], H[1], whichPass, 
                      )

Note that the last dimension of the arrays is inferred and thus not specified.