  module fft6_module

    use vtype_module

    implicit none

    private

    private :: fft6___fft_base2_core, fft6___fft_basearb_core

    public  :: fft6_fft_init, &
         fft6_fft, fft6_fftreal, fft6_fft1, fft6_fftreal1, &
         fft6_m2i, fft6_i2m, &
         fft6_fft_prototype, fft6_fft_base2


    !**************************************************************************

  contains

    !**************************************************************************
    !**************************************************************************
    ! mm      : Number of data set
    ! nn      : Data size, nn must be 2^x where x is integer
    ! ita     : Pre-calculated integer data for FFT
    ! dta     : Pre-calculated double-precision data for FFT
    ! datar   : Real part of data (input/output)
    ! datai   : Imaginary part of data (input/output)
    ! isign   : isign =  1
    !         : isign = -1     Inverse transformation
    !**************************************************************************

    subroutine fft6_fft_init( nn, ita, dta )

      use const_module, only : pix2

      integer(i4b), intent(in ) :: nn
      integer(i4b), intent(out) :: ita( 5 )
      real(dp)    , intent(out) :: dta( nn, 2 )


      !
      ! local variables
      !
      integer(i4b) :: basenum, np, i, j, l, tfindex
      real(dp)     :: theta, sinth, costh

      !
      ! bn   ( i ) :: base number for each ith-base FFT
      ! nloop( i ) :: number of loop for each ith-base FFT
      !
      integer(i4b), parameter :: nbn = 4
      integer(i4b)            :: bn( nbn ), nloop( nbn )

      !
      ! tf         :: twiddle factors
      ! tf( :, 1 ) :: sine   values for ith-base FFT
      ! tf( :, 2 ) :: cosine values for ith-base FFT
      !            :: i = 1 : 2-base FFT
      !            :: i = 2 : 3-base FFT
      !            :: i = 3 : 5-base FFT
      !
      real(dp)                :: tf( nn, 2 )



      bn( 1 ) =  2
      bn( 2 ) =  3
      bn( 3 ) =  5

      np = nn
      !
      do j = 1, nbn-1
         basenum = bn( j )
         l = 0
         find_base : do
            if( mod( np, basenum ) .ne. 0 ) exit find_base
            np = np / basenum
            l = l + 1
         end do find_base
         nloop( j ) = l
      end do
      !
      if( np .gt. 1 ) then
         write( 6, * ) 'WARNING:Unable to factorize the number ', np, &
              'by 2, 3, and 5.'
         bn   ( nbn ) = np
         nloop( nbn ) = 1
      else
         bn   ( nbn ) = 0
         nloop( nbn ) = 0
      end if



!!$      !
!!$      ! calculation of twiddle factor for 2-base FFT
!!$      !
!!$      tf( :, : ) = 0.0d0
!!$      tfindex    = 0
!!$      np         = nn
!!$      !
!!$      basenum = bn( 1 )
!!$      do l = 1, nloop( 1 )
!!$         np = np / basenum
!!$
!!$         do iph = 0, np - 1
!!$            tfindex = tfindex + 1
!!$            if( tfindex .gt. size( tf, 1 ) ) &
!!$                 stop 'Size of array "tf" is not enough.'
!!$            theta = pix2 * dble( iph ) / ( np + np )
!!$            sinth = sin( theta )
!!$            costh = cos( theta )
!!$            tf( tfindex, 1 ) = sinth
!!$            tf( tfindex, 2 ) = costh
!!$         end do
!!$
!!$      end do

      !
      ! calculation of twiddle factor
      !
      do i = 0, nn-1
         theta = pix2 * dble( i ) / nn
         sinth = sin( theta )
         costh = cos( theta )
         tf( i+1, 1 ) = sinth
         tf( i+1, 2 ) = costh
      end do



      do j = 1, nbn
         ita( j ) = nloop( j )
      end do
      ita( nbn+1 ) = bn( nbn )

      dta( :, : ) = tf( :, : )


    end subroutine fft6_fft_init

    !**************************************************************************
    !**************************************************************************
    ! The unit of transformed output would be the same as that of input data.
    !**************************************************************************
    ! Spectra are divided by nn (data size) in transforming from the physical
    ! space to spectral space.
    ! That is not the case in the inverse-transformation.
    !
    ! Using this program, the power spectrum density (energy per unit
    ! frequency), p_i, can be obtained using output of this routine, datar and
    ! datai, as follows:
    !   p_i = nn * dt * ( datar_i**2 + datai_i**2 ),
    ! where dt is the time increment.
    ! If you consider about the power spectrum density of space-series, the
    ! energy density should be defined as follows:
    !   p_i = nn * dx * ( datar_i**2 + datai_i**2 ) / ( 2 pi ),
    ! where dx is the distance between two points.
    ! This is because we normally use the unit of "angular wavenumber"
    ! expressing the wavenumber.
    !**************************************************************************
    ! Array exchange routine was deleted at 2001/03/21.
    ! The wavenumber corresponding to array(i) (index, i) is ordered
    ! as follows,
    ! 0,1, ... , nn/2-2, nn/2-1 nn/2,-(nn/2-1),-(nn/2-2), ..., -1 cycles.
    !**************************************************************************
    ! The definition of fourier transform is 
    ! <A> = int A(x) exp(-ikx) dk.
    !**************************************************************************
    ! mm      : Number of data set
    ! nn      : Data size, nn must be 2^x where x is integer
    ! ita     : Pre-calculated integer data for FFT
    ! dta     : Pre-calculated double-precision data for FFT
    ! datar   : Real part of data (input/output)
    ! datai   : Imaginary part of data (input/output)
    ! isign   : isign =  1
    !         : isign = -1     Inverse transformation
    !**************************************************************************

    subroutine fft6_fft( mm, nn, ita, dta, datar, datai, isign )

      integer(i4b), intent(in   ) :: mm, nn, ita( 5 )
      real(dp)    , intent(in   ) :: dta( nn, 2 )
      real(dp)    , intent(inout) :: datar( mm, nn ), datai( mm, nn )
      integer(i4b), intent(in   ) :: isign


      !
      ! local variables
      !
      !
      ! bn( i )       :: base number for each ith-base FFT
      ! nloop( i )    :: number of loop for each ith-base FFT
      !
      integer(i4b), parameter :: nbn = 4
      integer(i4b)            :: bn( nbn ), nloop( nbn )
      real(dp)                :: tf( nn, 2 )

      integer(i4b) :: np, j


      bn( 1   ) =  2
      bn( 2   ) =  3
      bn( 3   ) =  5
      bn( nbn ) = ita( nbn+1 )

      do j = 1, nbn
         nloop( j ) = ita( j )
      end do

      tf( :, : ) = dta( :, : )

      np = nn
      call fft6___fft_base2_core  ( mm, nn, datar, datai, np, tf, nloop( 1 ), isign )

      call fft6___fft_basearb_core( mm, nn, datar, datai, np, nbn, bn, nloop, isign )


      if( isign .eq. 1 ) then
         datar( :, : ) = datar( :, : ) / nn
         datai( :, : ) = datai( :, : ) / nn
      end if


    end subroutine fft6_fft

    !**************************************************************************

    subroutine fft6___fft_base2_core( mm, nn, datar, datai, np, tf, nloop, isign )

      integer(i4b), intent(in   ) :: mm, nn
      real(dp)    , intent(inout) :: datar( mm, nn ), datai( mm, nn )
      integer(i4b), intent(inout) :: np
      real(dp)    , intent(in   ) :: tf( nn, 2 )
      integer(i4b), intent(in   ) :: nloop
      integer(i4b), intent(in   ) :: isign


      !
      ! local variables
      !
      integer(i4b) :: basenum = 2
      real(dp)     :: data2r( mm, nn ), data2i( mm, nn )
      integer(i4b) :: nq, nr, iph, iqh, irh, iq
      real(dp)     :: sinth, costh
      integer(i4b) :: i0, i1, k0, k1
      integer(i4b) :: tfi, tfi_interval
      integer(i4b) :: l



      tfi_interval = 1

      do l = 0, nloop-1

         np  = np / basenum
         nq  = basenum
         nr  = nn / ( np * nq )

         tfi = 1

         do iph = 0, np - 1
            sinth = tf( tfi, 1 ) * ( -isign )
            costh = tf( tfi, 2 )

            do irh = 0, nr - 1

               iq = 0
               k0 = np * nr * iq + nr * iph + irh + 1
               iq = 1
               k1 = np * nr * iq + nr * iph + irh + 1

               iqh = 0
               i0 = nq * nr * iph + nr * iqh + irh + 1
               iqh = 1
               i1 = nq * nr * iph + nr * iqh + irh + 1

               data2r( :, i0 ) = datar( :, k0 ) + datar( :, k1 )
               data2i( :, i0 ) = datai( :, k0 ) + datai( :, k1 )

               data2r( :, i1 ) = ( datar( :, k0 ) - datar( :, k1 ) ) * costh &
                    - ( datai( :, k0 ) - datai( :, k1 ) ) * sinth
               data2i( :, i1 ) = ( datar( :, k0 ) - datar( :, k1 ) ) * sinth &
                    + ( datai( :, k0 ) - datai( :, k1 ) ) * costh

            end do

            tfi = tfi + tfi_interval
         end do

         datar( :, : ) = data2r( :, : )
         datai( :, : ) = data2i( :, : )

         tfi_interval = tfi_interval * basenum
      end do


    end subroutine fft6___fft_base2_core

    !**************************************************************************

    subroutine fft6___fft_basearb_core( mm, nn, datar, datai, np, nbn, bn, nloop, isign )

      use const_module, only : pix2

      integer(i4b), intent(in   ) :: mm, nn
      real(dp)    , intent(inout) :: datar( mm, nn ), datai( mm, nn )
      integer(i4b), intent(inout) :: np
      integer(i4b), intent(in   ) :: nbn
      integer(i4b), intent(in   ) :: bn( nbn ), nloop( nbn )
      integer(i4b), intent(in   ) :: isign


      !
      ! local variables
      !
      integer(i4b) :: basenum
      real(dp)     :: data2r( mm, nn ), data2i( mm, nn )
      integer(i4b) :: nq, nr, iph, iqh, irh, iq
      real(dp)     :: theta, sinth, costh
      integer(i4b) :: i, j, k, l


      loop_for_base : do j = 1+1, nbn

         do l = 0, nloop( j )-1
            basenum = bn( j )

            np = np / basenum
            nq = basenum
            nr = nn / ( np * nq )

            data2r( :, : ) = 0.0d0
            data2i( :, : ) = 0.0d0

            do iph = 0, np - 1

               do irh = 0, nr - 1

                  do iqh = 0, nq - 1

                     i = nq * nr * iph + nr * iqh + irh + 1

                     do iq = 0, nq - 1
                        theta = pix2 * dble( iqh ) / nq &
                             * ( dble( iph ) / np + iq )
                        sinth = sin( theta ) * ( -isign )
                        costh = cos( theta )

                        k = np * nr * iq + nr * iph + irh + 1

                        data2r( :, i ) = data2r( :, i ) &
                             + datar( :, k ) * costh - datai( :, k ) * sinth
                        data2i( :, i ) = data2i( :, i ) &
                             + datar( :, k ) * sinth + datai( :, k ) * costh

                     end do

                  end do

               end do
            end do

            datar( :, : ) = data2r( :, : )
            datai( :, : ) = data2i( :, : )

         end do

      end do loop_for_base


    end subroutine fft6___fft_basearb_core

    !**************************************************************************

    subroutine fft6_fftreal( mm, nn, ita, dta, data, fcr, fci, isign )

      integer(i4b), intent(in   ) :: mm, nn, ita( 5 )
      real(dp)    , intent(in   ) :: dta( nn, 2 )
      real(dp)    , intent(inout) :: data( mm, nn )
      real(dp)    , intent(inout) :: fcr( mm, nn ), fci( mm, nn )
      integer(i4b), intent(in   ) :: isign


      !
      ! local variables
      !
      real(dp)     :: &
!           data1r( max( mm/2, 1 ), nn ), data2r( max( mm/2, 1 ), nn ), &
           fc1r  ( max( mm/2, 1 ), nn ), fc1i  ( max( mm/2, 1 ), nn ), &
           fc2r  ( max( mm/2, 1 ), nn ), fc2i  ( max( mm/2, 1 ), nn )
      real(dp)     :: &
!           data1d( 1, nn ), &
           fc1dr( 1, nn ), fc1di( 1, nn )

      real(dp)     :: rep, rem, aip, aim

      integer(i4b) :: odd_sw

      integer(i4b) :: j, l, m, n, nn2


      odd_sw = mod( mm, 2 )


      if( isign .eq. 1 ) then

         if( mm .ge. 2 ) then

            !
            ! packing two real data into one complex data
            !
!---------------------------------------------------------------------
!            do n = 1, nn
!               do m = 1, mm / 2
!                  data1r( m, n ) = data( m     , n )
!                  data2r( m, n ) = data( m+mm/2, n )
!               enddo
!            enddo
!            do n = 1, nn
!               do m = 1, mm / 2
!                  fc1r( :, : ) = data1r( :, : )
!                  fc1i( :, : ) = data2r( :, : )
!               enddo
!            enddo
!---------------------------------------------------------------------
            do n = 1, nn
               do m = 1, mm / 2
                  fc1r( m, n ) = data( m     , n )
                  fc1i( m, n ) = data( m+mm/2, n )
               end do
            end do
!---------------------------------------------------------------------

            call fft6_fft( mm/2, nn, ita, dta, fc1r, fc1i, isign )

            !
            ! decomposition of one complex data to two complex data
            !
            nn2 = 1 + 1 + nn
            !
            do l = 1, mm/2
               fc2r( l, 1 ) = fc1i( l, 1 )
               fc1i( l, 1 ) = 0.0d0
               fc2i( l, 1 ) = 0.0d0
            end do
            do j = 2, nn/2+1
               do l = 1, mm/2
                  rep = 0.5d0 * ( fc1r( l, j ) + fc1r( l, nn2 - j ) )
                  rem = 0.5d0 * ( fc1r( l, j ) - fc1r( l, nn2 - j ) )
                  aip = 0.5d0 * ( fc1i( l, j ) + fc1i( l, nn2 - j ) )
                  aim = 0.5d0 * ( fc1i( l, j ) - fc1i( l, nn2 - j ) )
                  fc1r( l, j ) = rep
                  fc1i( l, j ) = aim
                  fc1r( l, nn2 - j ) = rep
                  fc1i( l, nn2 - j ) = -aim
                  fc2r( l, j ) = aip
                  fc2i( l, j ) = -rem
                  fc2r( l, nn2 - j ) = aip
                  fc2i( l, nn2 - j ) = rem
               end do
            end do

            do n = 1, nn
               do m = 1, mm / 2
                  fcr( m     , n ) = fc1r( m, n )
                  fcr( m+mm/2, n ) = fc2r( m, n )
                  fci( m     , n ) = fc1i( m, n )
                  fci( m+mm/2, n ) = fc2i( m, n )
               end do
            end do

         end if

         if( odd_sw .eq. 1 ) then
            do n = 1, nn
               fc1dr( 1, n ) = data( mm, n )
               fc1di( 1, n ) = 0.0d0
            end do

            call fft6_fft( 1, nn, ita, dta, fc1dr, fc1di, isign )

            do n = 1, nn
               fcr( mm, n ) = fc1dr( 1, n )
               fci( mm, n ) = fc1di( 1, n )
            end do
         end if


      else if( isign .eq. -1 ) then

         if( mm .ge. 2 ) then

            !
            ! add the second complex array multiplied with the imaginary unit
            ! to the first complex array
            !
!---------------------------------------------------------------------
!            do n = 1, nn
!               do m = 1, mm/ 2
!                  fc1r( m, n ) = fcr( m     , n )
!                  fc2r( m, n ) = fcr( m+mm/2, n )
!                  fc1i( m, n ) = fci( m     , n )
!                  fc2i( m, n ) = fci( m+mm/2, n )
!               enddo
!            enddo
!            fc1r( :, : ) = fc1r( :, : ) - fc2i( :, : )
!            fc1i( :, : ) = fc1i( :, : ) + fc2r( :, : )
!---------------------------------------------------------------------
            do n = 1, nn
               do m = 1, mm/ 2
                  fc1r( m, n ) = fcr( m, n ) - fci( m+mm/2, n )
                  fc1i( m, n ) = fci( m, n ) + fcr( m+mm/2, n )
               end do
            end do
!---------------------------------------------------------------------

            call fft6_fft( mm/2, nn, ita, dta, fc1r, fc1i, isign )

!---------------------------------------------------------------------
!            do n = 1, nn
!               do m = 1, mm/ 2
!                  data1r( m, n ) = fc1r( m, n )
!                  data2r( m, n ) = fc1i( m, n )
!               enddo
!            enddo
!            do n = 1, nn
!               do m = 1, mm/ 2
!                  data( m     , n ) = fc1r( m, n )
!                  data( m+mm/2, n ) = fc1i( m, n )
!               enddo
!            enddo
!---------------------------------------------------------------------
            do n = 1, nn
               do m = 1, mm/ 2
                  data( m     , n ) = fc1r( m, n )
                  data( m+mm/2, n ) = fc1i( m, n )
               end do
            end do
!---------------------------------------------------------------------

         end if

         if( odd_sw .eq. 1 ) then
            do n = 1, nn
               fc1dr( 1, n ) = fcr( mm, n )
               fc1di( 1, n ) = fci( mm, n )
            end do

            call fft6_fft( 1, nn, ita, dta, fc1dr, fc1di, isign )

            do n = 1, nn
               data( mm, n ) = fc1dr( 1, n )
            end do
         end if


      else
         write( 6, * ) 'STOP: Unexpected value isign in fftreal.'
         stop
      end if


    end subroutine fft6_fftreal

    !**************************************************************************

    subroutine fft6_fft1( nn, ita, dta, datar, datai, isign )

      integer(i4b), intent(in   ) :: nn, ita( 5 )
      real(dp)    , intent(in   ) :: dta( nn, 2 )
      real(dp)    , intent(inout) :: datar( nn ), datai( nn )
      integer(i4b), intent(in   ) :: isign


      !
      ! local variables
      !
      real(dp)     :: data2dr( 1, nn ), data2di( 1, nn )
      integer(i4b) :: n


      do n = 1, nn
         data2dr( 1, n ) = datar( n )
         data2di( 1, n ) = datai( n )
      end do

      call fft6_fft( 1, nn, ita, dta, data2dr, data2di, isign )

      do n = 1, nn
         datar( n ) = data2dr( 1, n )
         datai( n ) = data2di( 1, n )
      end do


    end subroutine fft6_fft1

    !**************************************************************************

    subroutine fft6_fftreal1( nn, ita, dta, data, fcr, fci, isign )

      integer(i4b), intent(in   ) :: nn, ita( 5 )
      real(dp)    , intent(in   ) :: dta( nn, 2 )
      real(dp)    , intent(inout) :: data( nn )
      real(dp)    , intent(inout) :: fcr( nn ), fci( nn )
      integer(i4b), intent(in   ) :: isign


      !
      ! local variables
      !
      real(dp)     :: data2d( 1, nn )
      real(dp)     :: fc2dr( 1, nn ), fc2di( 1, nn )
      integer(i4b) :: n


      select case( isign )
      case ( 1 )
         do n = 1, nn
            data2d( 1, n ) = data( n )
         end do
      case ( -1 )
         do n = 1, nn
            fc2dr ( 1, n ) = fcr ( n )
            fc2di ( 1, n ) = fci ( n )
         end do
      case default
         stop 'Unexpected isign in fft6_fftreal1.'
      end select

      call fft6_fftreal( 1, nn, ita, dta, data2d, fc2dr, fc2di, isign )

      select case( isign )
      case ( 1 )
         do n = 1, nn
            fcr( n ) = fc2dr( 1, n )
            fci( n ) = fc2di( 1, n )
         end do
      case ( -1 )
         do n = 1, nn
            data( n ) = data2d( 1, n )
         end do
      case default
         stop 'Unexpected isign in fft6_fftreal1.'
      end select


    end subroutine fft6_fftreal1

    !**************************************************************************
    ! function fft6_m2i
    !**************************************************************************
    ! nn    : Data size
    ! m     : Wavenumber
    !**************************************************************************

    function fft6_m2i( nn, m ) result( ires )

      integer(i4b)             :: ires

      integer(i4b), intent(in) :: nn, m


!      if( m .ge. 0 ) then
!         fft5_m2i = m + 1
!      else
!         fft5_m2i = nn + ( m + 1 )
!      endif


      ires = m + 1
      if( m .lt. 0 ) ires = nn + ires


!      fft5_m2i = ( min( m, 0 ) / m ) * nn + m + 1


    end function fft6_m2i

    !**************************************************************************
    ! function fft6_i2m
    !**************************************************************************
    ! nn    : Data size
    ! i     : Index for array including Fourier coefficients
    !**************************************************************************

    function fft6_i2m( nn, i ) result( ires )

      integer(i4b)             :: ires

      integer(i4b), intent(in) :: nn, i


!      if( i .le. nn/2+1 ) then
!         fft5_i2m = i - 1
!      else
!         fft5_i2m = i - 1 - nn
!      endif


      ires = i - 1
      if( i .gt. nn/2+1 ) ires = ires - nn


    end function fft6_i2m

    !**************************************************************************

    subroutine fft6_fft_prototype( mm, nn, datar, datai, isign )

      use const_module, only : pix2

      integer(i4b), intent(in   ) :: mm, nn
      real(dp)    , intent(inout) :: datar( mm, nn ), datai( mm, nn )
      integer(i4b), intent(in   ) :: isign


      !
      ! local variables
      !
      real(dp)     :: data2r( mm, nn ), data2i( mm, nn )
      integer(i4b) :: basenum
      integer(i4b) :: np, nq, nr, iph, iqh, irh, iq
      real(dp)     :: theta, sinth, costh
      integer(i4b) :: i, k


      np = nn

      outer_loop : do
         if( np .eq. 1 ) exit outer_loop

         if( mod( np, 2 ) .eq. 0 ) then
            basenum = 2
         else if( mod( np, 3 ) .eq. 0 ) then
            basenum = 3
         else if( mod( np, 5 ) .eq. 0 ) then
            basenum = 5
         else
            write( 6, * ) 'WARNING:Unable to factorize the number ', np, &
                 'by 2, 3, and 5.'
            basenum = np
         end if

         np = np / basenum
         nq = basenum
         nr = nn / ( np * nq )

         data2r( :, : ) = 0.0d0
         data2i( :, : ) = 0.0d0

         do irh = 0, nr - 1
            do iph = 0, np - 1

               do iqh = 0, nq - 1

                  i = nq * nr * iph + nr * iqh + irh + 1

                  do iq = 0, nq - 1
                     theta = pix2 * dble( iqh ) / nq * ( dble( iph ) / np + iq )
                     sinth = sin( theta ) * ( -isign )
                     costh = cos( theta )

                     k = np * nr * iq + nr * iph + irh + 1

                     data2r( :, i ) = data2r( :, i ) &
                          + datar( :, k ) * costh - datai( :, k ) * sinth
                     data2i( :, i ) = data2i( :, i ) &
                          + datar( :, k ) * sinth + datai( :, k ) * costh

                  end do

               end do

            end do
         end do

         datar( :, : ) = data2r( :, : )
         datai( :, : ) = data2i( :, : )

      end do outer_loop

      if( isign .eq. 1 ) then
         datar( :, : ) = datar( :, : ) / nn
         datai( :, : ) = datai( :, : ) / nn
      end if


    end subroutine fft6_fft_prototype

    !**************************************************************************

    subroutine fft6_fft_base2( mm, nn, datar, datai, isign )

      use const_module, only : pix2

      integer(i4b), intent(in   ) :: mm, nn
      real(dp)    , intent(inout) :: datar( mm, nn ), datai( mm, nn )
      integer(i4b), intent(in   ) :: isign


      !
      ! local variables
      !
      real(dp)     :: data2r( mm, nn ), data2i( mm, nn )
      integer(i4b) :: basenum
      integer(i4b) :: np, nq, nr, iph, iqh, irh, iq
      real(dp)     :: theta, sinth, costh
      integer(i4b) :: i0, i1, k0, k1
      integer(i4b) :: tmpnn, loopn, l


      basenum = 2


      tmpnn = nn
      loopn = 0
      find_power : do
         if( tmpnn .eq. 1 ) exit find_power
         tmpnn = tmpnn / basenum
         loopn = loopn + 1
      end do find_power


      np = nn

      do l = 0, loopn-1

         np = np / basenum
         nq = basenum
         nr = nn / ( np * nq )

         do irh = 0, nr - 1
            do iph = 0, np - 1

               theta = pix2 * dble( iph ) / ( np + np )
               sinth = sin( theta ) * ( -isign )
               costh = cos( theta )

               iq = 0
               k0 = np * nr * iq + nr * iph + irh + 1
               iq = 1
               k1 = np * nr * iq + nr * iph + irh + 1

               iqh = 0
               i0 = nq * nr * iph + nr * iqh + irh + 1
               iqh = 1
               i1 = nq * nr * iph + nr * iqh + irh + 1

               data2r( :, i0 ) = datar( :, k0 ) + datar( :, k1 )
               data2i( :, i0 ) = datai( :, k0 ) + datai( :, k1 )

               data2r( :, i1 ) = ( datar( :, k0 ) - datar( :, k1 ) ) * costh &
                    - ( datai( :, k0 ) - datai( :, k1 ) ) * sinth
               data2i( :, i1 ) = ( datar( :, k0 ) - datar( :, k1 ) ) * sinth &
                    + ( datai( :, k0 ) - datai( :, k1 ) ) * costh

            end do
         end do

         datar( :, : ) = data2r( :, : )
         datai( :, : ) = data2i( :, : )

      end do

      if( isign .eq. 1 ) then
         datar( :, : ) = datar( :, : ) / nn
         datai( :, : ) = datai( :, : ) / nn
      end if


    end subroutine fft6_fft_base2

    !**************************************************************************

  end module fft6_module
