InvertFourierModes Subroutine

    subroutine InvertFourierModes(Nk,Vk,V)
      ! Computes velocities in physical space from fourier modes.
      use, intrinsic :: iso_c_binding
      implicit none
      integer,      intent(in)    :: Nk(3)                                     !! Size of fourier modes
      complex(wp),  intent(in)    :: Vk(:,:,:,:)                               !! Fourier modes
      type(eulerian_obj_r),        &
                    intent(inout) :: V(3)
      ! Work variables
      type(C_PTR)             :: plan       = C_NULL_PTR
      type(C_PTR)             :: p_in       = C_NULL_PTR
      type(C_PTR)             :: p_out      = C_NULL_PTR
      complex(C_DOUBLE_COMPLEX),  &
                      pointer :: in(:,:,:)  => null()
      real(C_DOUBLE), pointer :: out(:,:,:) => null()
      integer                 :: mi,mj,mk
      integer                 :: i,j,k,dir

      include 'fftw3.f03'

      ! Allocate memory using FFTW to ensure SIMD allignment
      ! Note that the output must include +2 padding cells in the
      ! third dimension. This is a requirement from FFTW.
      p_in  = fftw_alloc_complex(int(Nk(1)*Nk(2)*Nk(3),         C_SIZE_T))
      p_out = fftw_alloc_real(   int(2*(Ng(1)/2+1)*Ng(2)*Ng(3), C_SIZE_T))

      ! Convert C pointers to Fortran arrays
      call c_f_pointer(p_in,  in,  [Nk(1),Nk(2),Nk(3)])
      call c_f_pointer(p_out, out, [Ng(1),Ng(2),Ng(3)])

      ! Create plan - dimensions must be inverted per FFTW.
      plan = fftw_plan_dft_c2r_3d(Ng(3), Ng(2), Ng(1), in, out, FFTW_MEASURE)


      do dir=1,3

        do mk=1,Nk(3)
          do mj=1,Nk(2)
            do mi=1,Nk(1)
               in(mi,mj,mk) = cmplx(Vk(dir,mi,mj,mk),kind=C_DOUBLE_COMPLEX)
            end do
          end do
        end do

        ! Execute plan (build inverse)
        call fftw_execute_dft_c2r(plan, in, out)

        ! Normalize (FFTW does not automatically normalize)
        do k=1,Ng(3)
          do j=1,Ng(2)
            do i=1,Ng(1)
              V(dir)%cell(i,j,k) = real(out(i,j,k),wp)/real(Ng(1)*Ng(2)*Ng(3),wp)
            end do
          end do
        end do

      end do

      ! Clean up
      call dfftw_destroy_plan(plan)
      call fftw_free(p_out)
      call fftw_free(p_in)

      return
    end subroutine InvertFourierModes