PROGRAM FFT_3D_2Decomp_MPI
      use mpi
      use, intrinsic :: iso_c_binding
      use decomp_2d
      use decomp_2d_fft
      implicit none
      integer, parameter :: L = 128
      integer, parameter :: M = 128
      integer, parameter :: N = 128
      real, parameter :: prec = 1.0e-7
      integer :: p_row, p_col
      integer :: nx, ny, nz
      integer :: nglob, nlocin, nlocout
      integer, dimension(2) :: dims
      integer, dimension(3) :: fft_start, fft_end, fft_size
      real :: diff, cdiff, ccdiff, pi, twopi
      real(mytype), allocatable, dimension(:) :: sinx, siny, sinz
      real(mytype), allocatable, dimension(:,:,:) :: in, out2
      complex(mytype), allocatable, dimension(:,:,:) :: out
      real(mytype) :: dummy
      integer :: ierror, i,j,k, numproc, mype
      integer, dimension(3) :: sizex, sizez

! ===== Initialize
       call MPI_INIT(ierror)
       call MPI_Comm_size(MPI_COMM_WORLD, numproc, ierror) 
       call MPI_COMM_RANK(MPI_COMM_WORLD, mype, ierror)

       dims = 0
       call MPI_Dims_create(numproc, 2, dims, ierror)
       if (dims(1) .gt. dims(2)) then
         dims(1) = dims(2)
         dims(2) = numproc/dims(1)
       endif
       p_row = dims(1)
       p_col = dims(2)

       if (mype .eq. 0) write(*,*) 'Initializing 2Decomp&FFT...' 

       call decomp_2d_init(L,M,N,p_row,p_col)
       call decomp_2d_fft_init
       call MPI_Barrier(MPI_COMM_WORLD, ierror)

       if (mype .eq. 0) write(*,*) 'Initialize variables...' 

       pi = atan(1.0)*4.0;
       twopi = 2.0*pi;

       do i =1, 3
          sizex(i) = xend(i) - xstart(i) + 1
          sizez(i) = zend(i) - zstart(i) + 1
       end do

       nglob = L * M * N
       nlocin = sizex(1) * sizex(2) * sizex(3)
       nlocout = sizez(1) * sizez(2) * sizez(3)

       if (mype .eq. 0) write(*,*) 'Allocating array...' 

       call decomp_2d_fft_get_size(fft_start,fft_end,fft_size)
       allocate (in(xstart(1):xend(1),xstart(2):xend(2),xstart(3):xend(3)))
       allocate (out(fft_start(1):fft_end(1), fft_start(2):fft_end(2), fft_start(3):fft_end(3)))
       allocate (out2(xstart(1):xend(1),xstart(2):xend(2),xstart(3):xend(3)))

       allocate (sinx(L))
       allocate (siny(M))
       allocate (sinz(N))

! ===== each processor gets its local portion of global data =====
       if (mype .eq. 0) write(*,*) 'Initializing array...' 
       do i =1, L
         sinx(i) = sin(twopi*real(i)/L)
       end do
       do j =1, M
         siny(j) = sin(twopi*real(j)/M)
       end do
       do k =1, N
         sinz(k) = sin(twopi*real(k)/N)
       end do
       do k=xstart(3),xend(3)
          do j=xstart(2),xend(2)
             dummy = sinz(k) * siny(j)
             do i=xstart(1),xend(1)
                in(i,j,k) = dummy * sinx(i)
             end do
          end do
       end do

! ===== 3D forward FFT =====
       if (mype .eq. 0) write(*,*) 'Forward FFT...' 
       call decomp_2d_fft_3d(in, out)
       call MPI_Barrier(MPI_COMM_WORLD, ierror)

! ===== Normalizzation =====
       if (mype .eq. 0) write(*,*) 'Normalizzation...' 
       do k=fft_start(3),fft_end(3)
          do j=fft_start(2),fft_end(2)
             do i=fft_start(1),fft_end(1)
                out(i,j,k) = out(i,j,k)/nglob
             end do
          end do
       end do
       call MPI_Barrier(MPI_COMM_WORLD, ierror)

! ===== 3D backward =====
       if (mype .eq. 0) write(*,*) 'Backward FFT...' 
       call decomp_2d_fft_3d(out, out2)

! ===== Compare =====
       call MPI_Barrier(MPI_COMM_WORLD, ierror)
       if (mype .eq. 0) write(*,*) 'Comparing...' 
       cdiff =0.
       do k=xstart(3),xend(3)
          do j=xstart(2),xend(2)
             do i=xstart(1),xend(1)
               diff = abs(in(i,j,k) - out2(i,j,k))
               if(cdiff .lt. diff) then
                 cdiff = diff
               endif
             end do
          end do
        end do
        call MPI_Reduce(cdiff, ccdiff, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD, ierror)
        if (mype .eq. 0) then
          write(*,*) 'Cdiff_max =', ccdiff
          if(cdiff .gt. (prec*nglob*0.25)) then
            write(*,*) 'Results are incorrect:'
          else
            write(*,*) 'Results are correct!'
          endif
        endif

! =========== Finalize ===============
       if (mype .eq. 0) write(*,*) 'Finalizzing...' 
        call decomp_2d_fft_finalize
        call decomp_2d_finalize
        deallocate(in, out, out2)
        call MPI_FINALIZE(ierror)

      end