! Fast Fourier/Cosine/Sine Transform
!     dimension   :three
!     data length :power of 2
!     decimation  :frequency
!     radix       :split-radix, row-column
!     data        :inplace
!     table       :use
! subroutines
!     cdft3d: Complex Discrete Fourier Transform
!     rdft3d: Real Discrete Fourier Transform
!     ddct3d: Discrete Cosine Transform
!     ddst3d: Discrete Sine Transform
! necessary package
!     fftsg.f  : 1D-FFT package
!
!
! -------- Complex DFT (Discrete Fourier Transform) --------
!     [definition]
!         <case1>
!             X(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               x(j1,j2,j3) * 
!                               exp(2*pi*i*j1*k1/n1) * 
!                               exp(2*pi*i*j2*k2/n2) * 
!                               exp(2*pi*i*j3*k3/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         <case2>
!             X(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               x(j1,j2,j3) * 
!                               exp(-2*pi*i*j1*k1/n1) * 
!                               exp(-2*pi*i*j2*k2/n2) * 
!                               exp(-2*pi*i*j3*k3/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call cdft3d(n1max, n2max, 2*n1, n2, n3, 1, a, t, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call cdft3d(n1max, n2max, 2*n1, n2, n3, -1, a, t, ip, w)
!     [parameters]
!         n1max  :row1 size of the 3D array (integer)
!         n2max  :row2 size of the 3D array (integer)
!         2*n1   :data length (integer)
!                 n1 >= 1, n1 = power of 2
!         n2     :data length (integer)
!                 n2 >= 1, n2 = power of 2
!         n3     :data length (integer)
!                 n3 >= 1, n3 = power of 2
!         a(0:2*n1-1,0:n2-1,0:n3-1)
!                :input/output data (real*8)
!                 input data
!                     a(2*j1,j2,j3) = Re(x(j1,j2,j3)), 
!                     a(2*j1+1,j2,j3) = Im(x(j1,j2,j3)), 
!                     0<=j1<n1, 0<=j2<n2, 0<=j3<n3
!                 output data
!                     a(2*k1,k2,k3) = Re(X(k1,k2,k3)), 
!                     a(2*k1+1,k2,k3) = Im(X(k1,k2,k3)), 
!                     0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         t(0:*) :work area (real*8)
!                 length of t >= max(8*n2, 8*n3)
!         ip(0:*):work area for bit reversal (integer)
!                 length of ip >= 2+sqrt(n)
!                 (n = max(n1, n2, n3))
!                 ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:*) :cos/sin table (real*8)
!                 length of w >= max(n1/2, n2/2, n3/2)
!                 w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of 
!             call cdft3d(n1max, n2max, 2*n1, n2, n3, -1, a, t, ip, w)
!         is 
!             call cdft3d(n1max, n2max, 2*n1, n2, n3, 1, a, t, ip, w)
!             do j3 = 0, n3 - 1
!                 do j2 = 0, n2 - 1
!                     do j1 = 0, 2 * n1 - 1
!                         a(j1,j2,j3) = a(j1,j2,j3) * (1.0d0/n1/n2/n3)
!                     end do
!                 end do
!             end do
!         .
!
!
! -------- Real DFT / Inverse of Real DFT --------
!     [definition]
!         <case1> RDFT
!             R(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               a(j1,j2,j3) * 
!                               cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + 
!                                   2*pi*j3*k3/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!             I(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               a(j1,j2,j3) * 
!                               sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + 
!                                   2*pi*j3*k3/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         <case2> IRDFT (excluding scale)
!             a(k1,k2,k3) = (1/2) * sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               (R(j1,j2,j3) * 
!                               cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + 
!                                   2*pi*j3*k3/n3) + 
!                               I(j1,j2,j3) * 
!                               sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2 + 
!                                   2*pi*j3*k3/n3)), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         (notes: R(mod(n1-k1,n1),mod(n2-k2,n2),mod(n3-k3,n3)) = R(k1,k2,k3), 
!                 I(mod(n1-k1,n1),mod(n2-k2,n2),mod(n3-k3,n3)) = -I(k1,k2,k3), 
!                 0<=k1<n1, 0<=k2<n2, 0<=k3<n3)
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call rdft3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call rdft3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!     [parameters]
!         n1max  :row1 size of the 3D array (integer)
!         n2max  :row2 size of the 3D array (integer)
!         n1     :data length (integer)
!                 n1 >= 2, n1 = power of 2
!         n2     :data length (integer)
!                 n2 >= 2, n2 = power of 2
!         n3     :data length (integer)
!                 n3 >= 2, n3 = power of 2
!         a(0:n1-1,0:n2-1,0:n3-1)
!                :input/output data (real*8)
!                 <case1>
!                     output data
!                         a(2*k1,k2,k3) = R(k1,k2,k3)
!                                 = R(n1-k1,mod(n2-k2,n2),mod(n3-k3,n3)), 
!                         a(2*k1+1,k2,k3) = I(k1,k2,k3)
!                                 = -I(n1-k1,mod(n2-k2,n2),mod(n3-k3,n3)), 
!                             0<k1<n1/2, 0<=k2<n2, 0<=k3<n3, 
!                         a(0,k2,k3) = R(0,k2,k3)
!                                    = R(0,n2-k2,mod(n3-k3,n3)), 
!                         a(1,k2,k3) = I(0,k2,k3)
!                                    = -I(0,n2-k2,mod(n3-k3,n3)), 
!                         a(1,n2-k2,k3) = R(n1/2,k2,mod(n3-k3,n3))
!                                       = R(n1/2,n2-k2,k3), 
!                         a(0,n2-k2,k3) = -I(n1/2,k2,mod(n3-k3,n3))
!                                       = I(n1/2,n2-k2,k3), 
!                             0<k2<n2/2, 0<=k3<n3, 
!                         a(0,0,k3) = R(0,0,k3)
!                                   = R(0,0,n3-k3), 
!                         a(1,0,k3) = I(0,0,k3)
!                                   = -I(0,0,n3-k3), 
!                         a(0,n2/2,k3) = R(0,n2/2,k3)
!                                      = R(0,n2/2,n3-k3), 
!                         a(1,n2/2,k3) = I(0,n2/2,k3)
!                                      = -I(0,n2/2,n3-k3), 
!                         a(1,0,n3-k3) = R(n1/2,0,k3)
!                                      = R(n1/2,0,n3-k3), 
!                         a(0,0,n3-k3) = -I(n1/2,0,k3)
!                                      = I(n1/2,0,n3-k3), 
!                         a(1,n2/2,n3-k3) = R(n1/2,n2/2,k3)
!                                         = R(n1/2,n2/2,n3-k3), 
!                         a(0,n2/2,n3-k3) = -I(n1/2,n2/2,k3)
!                                         = I(n1/2,n2/2,n3-k3), 
!                             0<k3<n3/2, 
!                         a(0,0,0) = R(0,0,0), 
!                         a(1,0,0) = R(n1/2,0,0), 
!                         a(0,0,n3/2) = R(0,0,n3/2), 
!                         a(1,0,n3/2) = R(n1/2,0,n3/2), 
!                         a(0,n2/2,0) = R(0,n2/2,0), 
!                         a(1,n2/2,0) = R(n1/2,n2/2,0), 
!                         a(0,n2/2,n3/2) = R(0,n2/2,n3/2), 
!                         a(1,n2/2,n3/2) = R(n1/2,n2/2,n3/2)
!                 <case2>
!                     input data
!                         a(2*j1,j2,j3) = R(j1,j2,j3)
!                                 = R(n1-j1,mod(n2-j2,n2),mod(n3-j3,n3)), 
!                         a(2*j1+1,j2,j3) = I(j1,j2,j3)
!                                 = -I(n1-j1,mod(n2-j2,n2),mod(n3-j3,n3)), 
!                             0<j1<n1/2, 0<=j2<n2, 0<=j3<n3, 
!                         a(0,j2,j3) = R(0,j2,j3)
!                                    = R(0,n2-j2,mod(n3-j3,n3)), 
!                         a(1,j2,j3) = I(0,j2,j3)
!                                    = -I(0,n2-j2,mod(n3-j3,n3)), 
!                         a(1,n2-j2,j3) = R(n1/2,j2,mod(n3-j3,n3))
!                                       = R(n1/2,n2-j2,j3), 
!                         a(0,n2-j2,j3) = -I(n1/2,j2,mod(n3-j3,n3))
!                                       = I(n1/2,n2-j2,j3), 
!                             0<j2<n2/2, 0<=j3<n3, 
!                         a(0,0,j3) = R(0,0,j3)
!                                   = R(0,0,n3-j3), 
!                         a(1,0,j3) = I(0,0,j3)
!                                   = -I(0,0,n3-j3), 
!                         a(0,n2/2,j3) = R(0,n2/2,j3)
!                                      = R(0,n2/2,n3-j3), 
!                         a(1,n2/2,j3) = I(0,n2/2,j3)
!                                      = -I(0,n2/2,n3-j3), 
!                         a(1,0,n3-j3) = R(n1/2,0,j3)
!                                      = R(n1/2,0,n3-j3), 
!                         a(0,0,n3-j3) = -I(n1/2,0,j3)
!                                      = I(n1/2,0,n3-j3), 
!                         a(1,n2/2,n3-j3) = R(n1/2,n2/2,j3)
!                                         = R(n1/2,n2/2,n3-j3), 
!                         a(0,n2/2,n3-j3) = -I(n1/2,n2/2,j3)
!                                         = I(n1/2,n2/2,n3-j3), 
!                             0<j3<n3/2, 
!                         a(0,0,0) = R(0,0,0), 
!                         a(1,0,0) = R(n1/2,0,0), 
!                         a(0,0,n3/2) = R(0,0,n3/2), 
!                         a(1,0,n3/2) = R(n1/2,0,n3/2), 
!                         a(0,n2/2,0) = R(0,n2/2,0), 
!                         a(1,n2/2,0) = R(n1/2,n2/2,0), 
!                         a(0,n2/2,n3/2) = R(0,n2/2,n3/2), 
!                         a(1,n2/2,n3/2) = R(n1/2,n2/2,n3/2)
!                 ---- output ordering ----
!                     call rdft3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!                     call rdft3dsort(n1max, n2max, n1, n2, n3, 1, a)
!                     ! stored data is a(0:n1-1,0:n2-1,0:n3+1):
!                     ! a(2*k1,k2,k3) = R(k1,k2,k3), 
!                     ! a(2*k1+1,k2,k3) = I(k1,k2,k3), 
!                     ! 0<=k1<=n1/2, 0<=k2<n2, 0<=k3<n3.
!                     ! the stored data is larger than the input data!
!                 ---- input ordering ----
!                     call rdft3dsort(n1max, n2max, n1, n2, n3, -1, a)
!                     call rdft3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!         t(0:*) :work area (real*8)
!                 length of t >= max(8*n2, 8*n3)
!         ip(0:*):work area for bit reversal (integer)
!                 length of ip >= 2+sqrt(n)
!                 (n = max(n1/2, n2, n3))
!                 ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:*) :cos/sin table (real*8)
!                 length of w >= max(n1/4, n2/2, n3/2) + n1/4
!                 w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of 
!             call rdft3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!         is 
!             call rdft3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!             do j3 = 0, n3 - 1
!                 do j2 = 0, n2 - 1
!                     do j1 = 0, n1 - 1
!                         a(j1,j2,j3) = a(j1,j2,j3) * (2.0d0/n1/n2/n3)
!                     end do
!                 end do
!             end do
!         .
!
!
! -------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
!     [definition]
!         <case1> IDCT (excluding scale)
!             C(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               a(j1,j2,j3) * 
!                               cos(pi*j1*(k1+1/2)/n1) * 
!                               cos(pi*j2*(k2+1/2)/n2) * 
!                               cos(pi*j3*(k3+1/2)/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         <case2> DCT
!             C(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               a(j1,j2,j3) * 
!                               cos(pi*(j1+1/2)*k1/n1) * 
!                               cos(pi*(j2+1/2)*k2/n2) * 
!                               cos(pi*(j3+1/2)*k3/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call ddct3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call ddct3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!     [parameters]
!         n1max  :row1 size of the 3D array (integer)
!         n2max  :row2 size of the 3D array (integer)
!         n1     :data length (integer)
!                 n1 >= 2, n1 = power of 2
!         n2     :data length (integer)
!                 n2 >= 2, n2 = power of 2
!         n3     :data length (integer)
!                 n3 >= 2, n3 = power of 2
!         a(0:n1-1,0:n2-1,0:n3-1)
!                :input/output data (real*8)
!                 output data
!                     a(k1,k2,k3) = C(k1,k2,k3), 
!                         0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         t(0:*) :work area (real*8)
!                 length of t >= max(4*n2, 4*n3)
!         ip(0:*):work area for bit reversal (integer)
!                 length of ip >= 2+sqrt(n)
!                 (n = max(n1/2, n2/2, n3/2))
!                 ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:*) :cos/sin table (real*8)
!                 length of w >= max(n1*3/2, n2*3/2, n3*3/2)
!                 w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of 
!             call ddct3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!         is 
!             do j3 = 0, n3 - 1
!                 do j2 = 0, n2 - 1
!                     a(0, j2, j3) = a(0, j2, j3) * 0.5d0
!                 end do
!                 do j1 = 0, n1 - 1
!                     a(j1, 0, j3) = a(j1, 0, j3) * 0.5d0
!                 end do
!             end do
!             do j2 = 0, n2 - 1
!                 do j1 = 0, n1 - 1
!                     a(j1, j2, 0) = a(j1, j2, 0) * 0.5d0
!                 end do
!             end do
!             call ddct3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!             do j3 = 0, n3 - 1
!                 do j2 = 0, n2 - 1
!                     do j1 = 0, n1 - 1
!                         a(j1,j2,j3) = a(j1,j2,j3) * (8.0d0/n1/n2/n3)
!                     end do
!                 end do
!             end do
!         .
!
!
! -------- DST (Discrete Sine Transform) / Inverse of DST --------
!     [definition]
!         <case1> IDST (excluding scale)
!             S(k1,k2,k3) = sum_j1=1^n1 sum_j2=1^n2 sum_j3=1^n3
!                               A(j1,j2,j3) * 
!                               sin(pi*j1*(k1+1/2)/n1) * 
!                               sin(pi*j2*(k2+1/2)/n2) * 
!                               sin(pi*j3*(k3+1/2)/n3), 
!                               0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!         <case2> DST
!             S(k1,k2,k3) = sum_j1=0^n1-1 sum_j2=0^n2-1 sum_j3=0^n3-1
!                               a(j1,j2,j3) * 
!                               sin(pi*(j1+1/2)*k1/n1) * 
!                               sin(pi*(j2+1/2)*k2/n2) * 
!                               sin(pi*(j3+1/2)*k3/n3), 
!                               0<k1<=n1, 0<k2<=n2, 0<k3<=n3
!     [usage]
!         <case1>
!             ip(0) = 0  ! first time only
!             call ddst3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!         <case2>
!             ip(0) = 0  ! first time only
!             call ddst3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!     [parameters]
!         n1max  :row1 size of the 3D array (integer)
!         n2max  :row2 size of the 3D array (integer)
!         n1     :data length (integer)
!                 n1 >= 2, n1 = power of 2
!         n2     :data length (integer)
!                 n2 >= 2, n2 = power of 2
!         n3     :data length (integer)
!                 n3 >= 2, n3 = power of 2
!         a(0:n1-1,0:n2-1,0:n3-1)
!                :input/output data (real*8)
!                 <case1>
!                     input data
!                         a(mod(j1,n1),mod(j2,n2),mod(j3,n3)) = A(j1,j2,j3), 
!                             0<j1<=n1, 0<j2<=n2, 0<j3<=n3
!                     output data
!                         a(k1,k2,k3) = S(k1,k2,k3), 
!                             0<=k1<n1, 0<=k2<n2, 0<=k3<n3
!                 <case2>
!                     output data
!                         a(mod(k1,n1),mod(k2,n2),mod(k3,n3)) = S(k1,k2,k3), 
!                             0<k1<=n1, 0<k2<=n2, 0<k3<=n3
!         t(0:*) :work area (real*8)
!                 length of t >= max(4*n2, 4*n3)
!         ip(0:*):work area for bit reversal (integer)
!                 length of ip >= 2+sqrt(n)
!                 (n = max(n1/2, n2/2, n3/2))
!                 ip(0),ip(1) are pointers of the cos/sin table.
!         w(0:*) :cos/sin table (real*8)
!                 length of w >= max(n1*3/2, n2*3/2, n3*3/2)
!                 w(),ip() are initialized if ip(0) = 0.
!     [remark]
!         Inverse of 
!             call ddst3d(n1max, n2max, n1, n2, n3, -1, a, t, ip, w)
!         is 
!             do j3 = 0, n3 - 1
!                 do j2 = 0, n2 - 1
!                     a(0, j2, j3) = a(0, j2, j3) * 0.5d0
!                 end do
!                 do j1 = 0, n1 - 1
!                     a(j1, 0, j3) = a(j1, 0, j3) * 0.5d0
!                 end do
!             end do
!             do j2 = 0, n2 - 1
!                 do j1 = 0, n1 - 1
!                     a(j1, j2, 0) = a(j1, j2, 0) * 0.5d0
!                 end do
!             end do
!             call ddst3d(n1max, n2max, n1, n2, n3, 1, a, t, ip, w)
!             do j3 = 0, n3 - 1
!                 do j2 = 0, n2 - 1
!                     do j1 = 0, n1 - 1
!                         a(j1,j2,j3) = a(j1,j2,j3) * (8.0d0/n1/n2/n3)
!                     end do
!                 end do
!             end do
!         .
!
!
      subroutine cdft3d(n1max, n2max, n1, n2, n3, isgn, a, 
     &    t, ip, w)
      integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), n
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      n = 2 * max(n2, n3)
      n = max(n, n1)
      if (n .gt. 4 * ip(0)) then
          call makewt(n / 4, ip, w)
      end if
      call xdft3da_sub(n1max, n2max, n1, n2, n3, 0, 
     &    isgn, a, t, ip, w)
      call cdft3db_sub(n1max, n2max, n1, n2, n3, 
     &    isgn, a, t, ip, w)
      end
!
      subroutine rdft3d(n1max, n2max, n1, n2, n3, isgn, a, 
     &    t, ip, w)
      integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), 
     &    n, nw, nc
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      n = 2 * max(n2, n3)
      n = max(n, n1)
      nw = ip(0)
      if (n .gt. 4 * nw) then
          nw = n / 4
          call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n1 .gt. 4 * nc) then
          nc = n1 / 4
          call makect(nc, ip, w(nw))
      end if
      if (isgn .lt. 0) then
          call rdft3d_sub(n1max, n2max, n1, n2, n3, isgn, a)
          call cdft3db_sub(n1max, n2max, n1, n2, n3, 
     &        isgn, a, t, ip, w)
          call xdft3da_sub(n1max, n2max, n1, n2, n3, 1, 
     &        isgn, a, t, ip, w)
      else
          call xdft3da_sub(n1max, n2max, n1, n2, n3, 1, 
     &        isgn, a, t, ip, w)
          call cdft3db_sub(n1max, n2max, n1, n2, n3, 
     &        isgn, a, t, ip, w)
          call rdft3d_sub(n1max, n2max, n1, n2, n3, isgn, a)
      end if
      end
!
      subroutine rdft3dsort(n1max, n2max, n1, n2, n3, isgn, a)
      integer n1max, n2max, n1, n2, n3, isgn, n2h, n3h, j, k
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), x, y
      n2h = n2 / 2
      n3h = n3 / 2
      if (isgn .lt. 0) then
          do k = 0, n3 - 1
              do j = n2h + 1, n2 - 1
                  a(0, j, k) = a(n1 + 1, j, k)
                  a(1, j, k) = a(n1, j, k)
              end do
          end do
          do k = n3h + 1, n3 - 1
              a(0, 0, k) = a(n1 + 1, 0, k)
              a(1, 0, k) = a(n1, 0, k)
              a(0, n2h, k) = a(n1 + 1, n2h, k)
              a(1, n2h, k) = a(n1, n2h, k)
          end do
          a(1, 0, 0) = a(n1, 0, 0)
          a(1, n2h, 0) = a(n1, n2h, 0)
          a(1, 0, n3h) = a(n1, 0, n3h)
          a(1, n2h, n3h) = a(n1, n2h, n3h)
      else
          do j = n2h + 1, n2 - 1
              y = a(0, j, 0)
              x = a(1, j, 0)
              a(n1, j, 0) = x
              a(n1 + 1, j, 0) = y
              a(n1, n2 - j, 0) = x
              a(n1 + 1, n2 - j, 0) = -y
              a(0, j, 0) = a(0, n2 - j, 0)
              a(1, j, 0) = -a(1, n2 - j, 0)
          end do
          do k = 1, n3 - 1
              do j = n2h + 1, n2 - 1
                  y = a(0, j, k)
                  x = a(1, j, k)
                  a(n1, j, k) = x
                  a(n1 + 1, j, k) = y
                  a(n1, n2 - j, n3 - k) = x
                  a(n1 + 1, n2 - j, n3 - k) = -y
                  a(0, j, k) = a(0, n2 - j, n3 - k)
                  a(1, j, k) = -a(1, n2 - j, n3 - k)
              end do
          end do
          do k = n3h + 1, n3 - 1
              y = a(0, 0, k)
              x = a(1, 0, k)
              a(n1, 0, k) = x
              a(n1 + 1, 0, k) = y
              a(n1, 0, n3 - k) = x
              a(n1 + 1, 0, n3 - k) = -y
              a(0, 0, k) = a(0, 0, n3 - k)
              a(1, 0, k) = -a(1, 0, n3 - k)
              y = a(0, n2h, k)
              x = a(1, n2h, k)
              a(n1, n2h, k) = x
              a(n1 + 1, n2h, k) = y
              a(n1, n2h, n3 - k) = x
              a(n1 + 1, n2h, n3 - k) = -y
              a(0, n2h, k) = a(0, n2h, n3 - k)
              a(1, n2h, k) = -a(1, n2h, n3 - k)
          end do
          a(n1, 0, 0) = a(1, 0, 0)
          a(n1 + 1, 0, 0) = 0
          a(1, 0, 0) = 0
          a(n1, n2h, 0) = a(1, n2h, 0)
          a(n1 + 1, n2h, 0) = 0
          a(1, n2h, 0) = 0
          a(n1, 0, n3h) = a(1, 0, n3h)
          a(n1 + 1, 0, n3h) = 0
          a(1, 0, n3h) = 0
          a(n1, n2h, n3h) = a(1, n2h, n3h)
          a(n1 + 1, n2h, n3h) = 0
          a(1, n2h, n3h) = 0
      end if
      end
!
      subroutine ddct3d(n1max, n2max, n1, n2, n3, isgn, a, 
     &    t, ip, w)
      integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), 
     &    n, nw, nc
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      n = max(n2, n3)
      n = max(n, n1)
      nw = ip(0)
      if (n .gt. 4 * nw) then
          nw = n / 4
          call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. nc) then
          nc = n
          call makect(nc, ip, w(nw))
      end if
      call ddxt3da_sub(n1max, n2max, n1, n2, n3, 0, 
     &    isgn, a, t, ip, w)
      call ddxt3db_sub(n1max, n2max, n1, n2, n3, 0, 
     &    isgn, a, t, ip, w)
      end
!
      subroutine ddst3d(n1max, n2max, n1, n2, n3, isgn, a, 
     &    t, ip, w)
      integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), 
     &    n, nw, nc
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      n = max(n2, n3)
      n = max(n, n1)
      nw = ip(0)
      if (n .gt. 4 * nw) then
          nw = n / 4
          call makewt(nw, ip, w)
      end if
      nc = ip(1)
      if (n .gt. nc) then
          nc = n
          call makect(nc, ip, w(nw))
      end if
      call ddxt3da_sub(n1max, n2max, n1, n2, n3, 1, 
     &    isgn, a, t, ip, w)
      call ddxt3db_sub(n1max, n2max, n1, n2, n3, 1, 
     &    isgn, a, t, ip, w)
      end
!
! -------- child routines --------
!
      subroutine xdft3da_sub(n1max, n2max, n1, n2, n3, icr, 
     &    isgn, a, t, ip, w)
      integer n1max, n2max, n1, n2, n3, icr, isgn, 
     &    ip(0 : *), i, j, k
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      do k = 0, n3 - 1
          if (icr .eq. 0) then
              do j = 0, n2 - 1
                  call cdft(n1, isgn, a(0, j, k), ip, w)
              end do
          else if (isgn .ge. 0) then
              do j = 0, n2 - 1
                  call rdft(n1, isgn, a(0, j, k), ip, w)
              end do
          end if
          if (n1 .gt. 4) then
              do i = 0, n1 - 8, 8
                  do j = 0, n2 - 1
                      t(2 * j) = a(i, j, k)
                      t(2 * j + 1) = a(i + 1, j, k)
                      t(2 * n2 + 2 * j) = a(i + 2, j, k)
                      t(2 * n2 + 2 * j + 1) = a(i + 3, j, k)
                      t(4 * n2 + 2 * j) = a(i + 4, j, k)
                      t(4 * n2 + 2 * j + 1) = a(i + 5, j, k)
                      t(6 * n2 + 2 * j) = a(i + 6, j, k)
                      t(6 * n2 + 2 * j + 1) = a(i + 7, j, k)
                  end do
                  call cdft(2 * n2, isgn, t, ip, w)
                  call cdft(2 * n2, isgn, t(2 * n2), ip, w)
                  call cdft(2 * n2, isgn, t(4 * n2), ip, w)
                  call cdft(2 * n2, isgn, t(6 * n2), ip, w)
                  do j = 0, n2 - 1
                      a(i, j, k) = t(2 * j)
                      a(i + 1, j, k) = t(2 * j + 1)
                      a(i + 2, j, k) = t(2 * n2 + 2 * j)
                      a(i + 3, j, k) = t(2 * n2 + 2 * j + 1)
                      a(i + 4, j, k) = t(4 * n2 + 2 * j)
                      a(i + 5, j, k) = t(4 * n2 + 2 * j + 1)
                      a(i + 6, j, k) = t(6 * n2 + 2 * j)
                      a(i + 7, j, k) = t(6 * n2 + 2 * j + 1)
                  end do
              end do
          else if (n1 .eq. 4) then
              do j = 0, n2 - 1
                  t(2 * j) = a(0, j, k)
                  t(2 * j + 1) = a(1, j, k)
                  t(2 * n2 + 2 * j) = a(2, j, k)
                  t(2 * n2 + 2 * j + 1) = a(3, j, k)
              end do
              call cdft(2 * n2, isgn, t, ip, w)
              call cdft(2 * n2, isgn, t(2 * n2), ip, w)
              do j = 0, n2 - 1
                  a(0, j, k) = t(2 * j)
                  a(1, j, k) = t(2 * j + 1)
                  a(2, j, k) = t(2 * n2 + 2 * j)
                  a(3, j, k) = t(2 * n2 + 2 * j + 1)
              end do
          else if (n1 .eq. 2) then
              do j = 0, n2 - 1
                  t(2 * j) = a(0, j, k)
                  t(2 * j + 1) = a(1, j, k)
              end do
              call cdft(2 * n2, isgn, t, ip, w)
              do j = 0, n2 - 1
                  a(0, j, k) = t(2 * j)
                  a(1, j, k) = t(2 * j + 1)
              end do
          end if
          if (icr .ne. 0 .and. isgn .lt. 0) then
              do j = 0, n2 - 1
                  call rdft(n1, isgn, a(0, j, k), ip, w)
              end do
          end if
      end do
      end
!
      subroutine cdft3db_sub(n1max, n2max, n1, n2, n3, 
     &    isgn, a, t, ip, w)
      integer n1max, n2max, n1, n2, n3, isgn, ip(0 : *), 
     &    i, j, k
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      if (n1 .gt. 4) then
          do j = 0, n2 - 1
              do i = 0, n1 - 8, 8
                  do k = 0, n3 - 1
                      t(2 * k) = a(i, j, k)
                      t(2 * k + 1) = a(i + 1, j, k)
                      t(2 * n3 + 2 * k) = a(i + 2, j, k)
                      t(2 * n3 + 2 * k + 1) = a(i + 3, j, k)
                      t(4 * n3 + 2 * k) = a(i + 4, j, k)
                      t(4 * n3 + 2 * k + 1) = a(i + 5, j, k)
                      t(6 * n3 + 2 * k) = a(i + 6, j, k)
                      t(6 * n3 + 2 * k + 1) = a(i + 7, j, k)
                  end do
                  call cdft(2 * n3, isgn, t, ip, w)
                  call cdft(2 * n3, isgn, t(2 * n3), ip, w)
                  call cdft(2 * n3, isgn, t(4 * n3), ip, w)
                  call cdft(2 * n3, isgn, t(6 * n3), ip, w)
                  do k = 0, n3 - 1
                      a(i, j, k) = t(2 * k)
                      a(i + 1, j, k) = t(2 * k + 1)
                      a(i + 2, j, k) = t(2 * n3 + 2 * k)
                      a(i + 3, j, k) = t(2 * n3 + 2 * k + 1)
                      a(i + 4, j, k) = t(4 * n3 + 2 * k)
                      a(i + 5, j, k) = t(4 * n3 + 2 * k + 1)
                      a(i + 6, j, k) = t(6 * n3 + 2 * k)
                      a(i + 7, j, k) = t(6 * n3 + 2 * k + 1)
                  end do
              end do
          end do
      else if (n1 .eq. 4) then
          do j = 0, n2 - 1
              do k = 0, n3 - 1
                  t(2 * k) = a(0, j, k)
                  t(2 * k + 1) = a(1, j, k)
                  t(2 * n3 + 2 * k) = a(2, j, k)
                  t(2 * n3 + 2 * k + 1) = a(3, j, k)
              end do
              call cdft(2 * n3, isgn, t, ip, w)
              call cdft(2 * n3, isgn, t(2 * n3), ip, w)
              do k = 0, n3 - 1
                  a(0, j, k) = t(2 * k)
                  a(1, j, k) = t(2 * k + 1)
                  a(2, j, k) = t(2 * n3 + 2 * k)
                  a(3, j, k) = t(2 * n3 + 2 * k + 1)
              end do
          end do
      else if (n1 .eq. 2) then
          do j = 0, n2 - 1
              do k = 0, n3 - 1
                  t(2 * k) = a(0, j, k)
                  t(2 * k + 1) = a(1, j, k)
              end do
              call cdft(2 * n3, isgn, t, ip, w)
              do k = 0, n3 - 1
                  a(0, j, k) = t(2 * k)
                  a(1, j, k) = t(2 * k + 1)
              end do
          end do
      end if
      end
!
      subroutine rdft3d_sub(n1max, n2max, n1, n2, n3, isgn, a)
      integer n1max, n2max, n1, n2, n3, isgn, 
     &    n2h, n3h, i, j, k, l
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), xi
      n2h = n2 / 2
      n3h = n3 / 2
      if (isgn .lt. 0) then
          do k = 1, n3h - 1
              l = n3 - k
              xi = a(0, 0, k) - a(0, 0, l)
              a(0, 0, k) = a(0, 0, k) + a(0, 0, l)
              a(0, 0, l) = xi
              xi = a(1, 0, l) - a(1, 0, k)
              a(1, 0, k) = a(1, 0, k) + a(1, 0, l)
              a(1, 0, l) = xi
              xi = a(0, n2h, k) - a(0, n2h, l)
              a(0, n2h, k) = a(0, n2h, k) + a(0, n2h, l)
              a(0, n2h, l) = xi
              xi = a(1, n2h, l) - a(1, n2h, k)
              a(1, n2h, k) = a(1, n2h, k) + a(1, n2h, l)
              a(1, n2h, l) = xi
              do i = 1, n2h - 1
                  j = n2 - i
                  xi = a(0, i, k) - a(0, j, l)
                  a(0, i, k) = a(0, i, k) + a(0, j, l)
                  a(0, j, l) = xi
                  xi = a(1, j, l) - a(1, i, k)
                  a(1, i, k) = a(1, i, k) + a(1, j, l)
                  a(1, j, l) = xi
                  xi = a(0, i, l) - a(0, j, k)
                  a(0, i, l) = a(0, i, l) + a(0, j, k)
                  a(0, j, k) = xi
                  xi = a(1, j, k) - a(1, i, l)
                  a(1, i, l) = a(1, i, l) + a(1, j, k)
                  a(1, j, k) = xi
              end do
          end do
          do i = 1, n2h - 1
              j = n2 - i
              xi = a(0, i, 0) - a(0, j, 0)
              a(0, i, 0) = a(0, i, 0) + a(0, j, 0)
              a(0, j, 0) = xi
              xi = a(1, j, 0) - a(1, i, 0)
              a(1, i, 0) = a(1, i, 0) + a(1, j, 0)
              a(1, j, 0) = xi
              xi = a(0, i, n3h) - a(0, j, n3h)
              a(0, i, n3h) = a(0, i, n3h) + a(0, j, n3h)
              a(0, j, n3h) = xi
              xi = a(1, j, n3h) - a(1, i, n3h)
              a(1, i, n3h) = a(1, i, n3h) + a(1, j, n3h)
              a(1, j, n3h) = xi
          end do
      else
          do k = 1, n3h - 1
              l = n3 - k
              a(0, 0, l) = 0.5d0 * (a(0, 0, k) - a(0, 0, l))
              a(0, 0, k) = a(0, 0, k) - a(0, 0, l)
              a(1, 0, l) = 0.5d0 * (a(1, 0, k) + a(1, 0, l))
              a(1, 0, k) = a(1, 0, k) - a(1, 0, l)
              a(0, n2h, l) = 0.5d0 * (a(0, n2h, k) - a(0, n2h, l))
              a(0, n2h, k) = a(0, n2h, k) - a(0, n2h, l)
              a(1, n2h, l) = 0.5d0 * (a(1, n2h, k) + a(1, n2h, l))
              a(1, n2h, k) = a(1, n2h, k) - a(1, n2h, l)
              do i = 1, n2h - 1
                  j = n2 - i
                  a(0, j, l) = 0.5d0 * (a(0, i, k) - a(0, j, l))
                  a(0, i, k) = a(0, i, k) - a(0, j, l)
                  a(1, j, l) = 0.5d0 * (a(1, i, k) + a(1, j, l))
                  a(1, i, k) = a(1, i, k) - a(1, j, l)
                  a(0, j, k) = 0.5d0 * (a(0, i, l) - a(0, j, k))
                  a(0, i, l) = a(0, i, l) - a(0, j, k)
                  a(1, j, k) = 0.5d0 * (a(1, i, l) + a(1, j, k))
                  a(1, i, l) = a(1, i, l) - a(1, j, k)
              end do
          end do
          do i = 1, n2h - 1
              j = n2 - i
              a(0, j, 0) = 0.5d0 * (a(0, i, 0) - a(0, j, 0))
              a(0, i, 0) = a(0, i, 0) - a(0, j, 0)
              a(1, j, 0) = 0.5d0 * (a(1, i, 0) + a(1, j, 0))
              a(1, i, 0) = a(1, i, 0) - a(1, j, 0)
              a(0, j, n3h) = 0.5d0 * (a(0, i, n3h) - a(0, j, n3h))
              a(0, i, n3h) = a(0, i, n3h) - a(0, j, n3h)
              a(1, j, n3h) = 0.5d0 * (a(1, i, n3h) + a(1, j, n3h))
              a(1, i, n3h) = a(1, i, n3h) - a(1, j, n3h)
          end do
      end if
      end
!
      subroutine ddxt3da_sub(n1max, n2max, n1, n2, n3, ics, 
     &    isgn, a, t, ip, w)
      integer n1max, n2max, n1, n2, n3, ics, isgn, 
     &    ip(0 : *), i, j, k
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      do k = 0, n3 - 1
          if (ics .eq. 0) then
              do j = 0, n2 - 1
                  call ddct(n1, isgn, a(0, j, k), ip, w)
              end do
          else
              do j = 0, n2 - 1
                  call ddst(n1, isgn, a(0, j, k), ip, w)
              end do
          end if
          if (n1 .gt. 2) then
              do i = 0, n1 - 4, 4
                  do j = 0, n2 - 1
                      t(j) = a(i, j, k)
                      t(n2 + j) = a(i + 1, j, k)
                      t(2 * n2 + j) = a(i + 2, j, k)
                      t(3 * n2 + j) = a(i + 3, j, k)
                  end do
                  if (ics .eq. 0) then
                      call ddct(n2, isgn, t, ip, w)
                      call ddct(n2, isgn, t(n2), ip, w)
                      call ddct(n2, isgn, t(2 * n2), ip, w)
                      call ddct(n2, isgn, t(3 * n2), ip, w)
                  else
                      call ddst(n2, isgn, t, ip, w)
                      call ddst(n2, isgn, t(n2), ip, w)
                      call ddst(n2, isgn, t(2 * n2), ip, w)
                      call ddst(n2, isgn, t(3 * n2), ip, w)
                  end if
                  do j = 0, n2 - 1
                      a(i, j, k) = t(j)
                      a(i + 1, j, k) = t(n2 + j)
                      a(i + 2, j, k) = t(2 * n2 + j)
                      a(i + 3, j, k) = t(3 * n2 + j)
                  end do
              end do
          else if (n1 .eq. 2) then
              do j = 0, n2 - 1
                  t(j) = a(0, j, k)
                  t(n2 + j) = a(1, j, k)
              end do
              if (ics .eq. 0) then
                  call ddct(n2, isgn, t, ip, w)
                  call ddct(n2, isgn, t(n2), ip, w)
              else
                  call ddst(n2, isgn, t, ip, w)
                  call ddst(n2, isgn, t(n2), ip, w)
              end if
              do j = 0, n2 - 1
                  a(0, j, k) = t(j)
                  a(1, j, k) = t(n2 + j)
              end do
          end if
      end do
      end
!
      subroutine ddxt3db_sub(n1max, n2max, n1, n2, n3, ics, 
     &    isgn, a, t, ip, w)
      integer n1max, n2max, n1, n2, n3, ics, isgn, 
     &    ip(0 : *), i, j, k
      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : n3 - 1), 
     &    t(0 : *), w(0 : *)
      if (n1 .gt. 2) then
          do j = 0, n2 - 1
              do i = 0, n1 - 4, 4
                  do k = 0, n3 - 1
                      t(k) = a(i, j, k)
                      t(n3 + k) = a(i + 1, j, k)
                      t(2 * n3 + k) = a(i + 2, j, k)
                      t(3 * n3 + k) = a(i + 3, j, k)
                  end do
                  if (ics .eq. 0) then
                      call ddct(n3, isgn, t, ip, w)
                      call ddct(n3, isgn, t(n3), ip, w)
                      call ddct(n3, isgn, t(2 * n3), ip, w)
                      call ddct(n3, isgn, t(3 * n3), ip, w)
                  else
                      call ddst(n3, isgn, t, ip, w)
                      call ddst(n3, isgn, t(n3), ip, w)
                      call ddst(n3, isgn, t(2 * n3), ip, w)
                      call ddst(n3, isgn, t(3 * n3), ip, w)
                  end if
                  do k = 0, n3 - 1
                      a(i, j, k) = t(k)
                      a(i + 1, j, k) = t(n3 + k)
                      a(i + 2, j, k) = t(2 * n3 + k)
                      a(i + 3, j, k) = t(3 * n3 + k)
                  end do
              end do
          end do
      else if (n1 .eq. 2) then
          do j = 0, n2 - 1
              do k = 0, n3 - 1
                  t(k) = a(0, j, k)
                  t(n3 + k) = a(1, j, k)
              end do
              if (ics .eq. 0) then
                  call ddct(n3, isgn, t, ip, w)
                  call ddct(n3, isgn, t(n3), ip, w)
              else
                  call ddst(n3, isgn, t, ip, w)
                  call ddst(n3, isgn, t(n3), ip, w)
              end if
              do k = 0, n3 - 1
                  a(0, j, k) = t(k)
                  a(1, j, k) = t(n3 + k)
              end do
          end do
      end if
      end
!