fft_obj_Transform Subroutine

    impure subroutine fft_obj_Transform(this,field)
      !> Computes the Fourier Transform of field and stores it back
      ! in the input field in the "halfcomplex" (r2r) format.
      ! E.g. if the FFT is to be carried out in the 1-dir, then
      ! field%cell(:,j,k) contains the fourier transform in the 1-dir
      ! stored as follows
      !    r_0, r_1, ..., r_(Ng(1)1/2), i_(Ng(1)+1/2-1),..., i_1
      !       r_m is the real part of mode m*2pi/L(1) and
      !       i_m is the imaginary part of mode m*2pi/L(1)
      ! for points at location y(j) (or ym(j), depending on the field
      ! staggering), z(k) (or zm(k)).
      implicit none
      class(fft_obj),       intent(inout) :: this                              !! A FFT object
      type(eulerian_obj_r), intent(inout) :: field                             !! Input field
      ! Work variables
      integer       :: n,m,dir
      integer       :: contig_dir_
      integer       :: Ng_T(3)
      integer       :: l_dims_T(3)
      integer       :: l_offsets_T(3)
      real(wp),      &
        allocatable :: array(:,:,:)
      real(wp),      &
        allocatable :: array2(:,:,:)

      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_f(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_f(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_f(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_Transform