fft_obj_Inverse Subroutine

    impure subroutine fft_obj_Inverse(this,field)
      !> Computes the Inverse of a Fourier Transform and stores it back
      ! in the input field. The FFT must be given in the halfcomplex
      ! format.
      implicit none
      class(fft_obj),       intent(inout) :: this                              !! A FFT object
      type(eulerian_obj_r), intent(inout) :: field                             !! Field containing FFT coefficient in halfcomplex format
      ! Work variables
      integer       :: n,m,dir
      integer       :: contig_dir_
      real(wp),      &
        allocatable :: array(:,:,:)
      real(wp),      &
        allocatable :: array2(:,:,:)
      integer       :: Ng_T(3)
      integer       :: l_dims_T(3)
      integer       :: l_offsets_T(3)

      do dir=1,3
        ! Check whether to take a Fourier transform in this direction
        if (this%fft_dir(dir)) then

          ! Check whether this direction is contiguous or decomposed
          if (this%parallel%np(dir).eq.1) then
            ! Direction is contiguous
            contig_dir_ = dir

          else
            ! Direction is decomposed
            contig_dir_ = this%contig_dir
          end if

          call this%Transpose(dir,contig_dir_,                   &
                   field%cell(this%block%lo(1):this%block%hi(1), &
                              this%block%lo(2):this%block%hi(2), &
                              this%block%lo(3):this%block%hi(3)),&
                   this%Ng,this%l_dims,this%l_offsets,           &
                   array, Ng_T, l_dims_T, l_offsets_T)

          select case (contig_dir_)
          case (1)

            do n=1,size(array,dim=3)
              do m=1,size(array,dim=2)
                this%buffers(dir)%in = real(array(:,m,n),C_DOUBLE)
                call fftw_execute_r2r(this%plan_b(dir),this%buffers(dir)%in,this%buffers(dir)%out)
                array(:,m,n) = real(this%buffers(dir)%out,wp)

              end do
            end do

          case (2)

            do n=1,size(array,dim=3)
              do m=1,size(array,dim=1)
                this%buffers(dir)%in = real(array(m,:,n),C_DOUBLE)
                call fftw_execute_r2r(this%plan_b(dir),this%buffers(dir)%in,this%buffers(dir)%out)
                array(m,:,n) = real(this%buffers(dir)%out,wp)
              end do
            end do

          case (3)

            do n=1,size(array,dim=2)
              do m=1,size(array,dim=1)
                this%buffers(dir)%in = real(array(m,n,:),C_DOUBLE)
                call fftw_execute_r2r(this%plan_b(dir),this%buffers(dir)%in,this%buffers(dir)%out)
                array(m,n,:) = real(this%buffers(dir)%out,wp)
              end do
            end do
          end select

          ! Transpose again
          call this%Transpose(dir,contig_dir_,            &
               array,Ng_T,l_dims_T,l_offsets_T,           &
               array2,this%Ng, this%l_dims, this%l_offsets)

          field%cell(this%block%lo(1):this%block%hi(1), &
                     this%block%lo(2):this%block%hi(2), &
                     this%block%lo(3):this%block%hi(3))=array2

        end if
      end do

      return
    end subroutine fft_obj_Inverse