subroutine da_solve_poissoneqn_fst(grid, xbx, del2b, b)

   !--------------------------------------------------------------------------
   !  Purpose: Solve Del**2 B = A for B with zero field boundary conditions.
   !
   !  Method:  1) Compute spectral del2b using double forward FST.
   !           2) Calculate spectral b.
   !           3) Reform gridpt. b using inverse double FST.
   !           Note no mean removed (would make b.ne.0 at boundaries) so
   !           arbitrary constant still present.
   !--------------------------------------------------------------------------

   implicit none

   type (domain),   intent(inout) :: grid
   type (xbx_type), intent(in)    :: xbx           ! Header & non-gridded vars.
   real,            intent(in)    :: del2b(ims:ime,jms:jme,kms:kme)  ! Del**2 B.
   real,            intent(out)   :: b(ims:ime,jms:jme,kms:kme)  ! B.
   
   integer           :: vector_inc    ! Increment between FST data.
   integer           :: vector_jump   ! Jump between start of vectors.
   integer           :: vector_size   ! Of form 2**p 3**q 5**r for FSTs.
   integer           :: num_vectors   ! Number of FSTs to perform.
   integer           :: work_area     ! Dimension for FST routine.
   integer           :: idim          ! Size of 1st dimension for FST.
   integer           :: jdim          ! Size of 2nd dimension for FST.

   integer           :: i, j, k, n, ij     ! loop counter

   real, allocatable :: work_1d(:)     ! FFT work array

   if (trace_use) call da_trace_entry("da_solve_poissoneqn_fst")

   !------------------------------------------------------------------------------
   ! [1.0] Initialise:
   !------------------------------------------------------------------------------

   ! Calculate work space needed.

   n = max(xbx%fft_ix*(grid%xp%jtex-grid%xp%jtsx+1), &
           xbx%fft_jy*(grid%xp%itey-grid%xp%itsy+1+xbx%pad_num))

   ! Allocate work arrays.
   allocate(work_1d(1:n))

   ! Copy del2b for transpose.

   grid%xp%v1z(its:ite,jts:jte,kts:kte) = del2b(its:ite,jts:jte,kts:kte)

   if (ite == ide) grid%xp%v1z(ite,jts:jte,kts:kte) = 0.0
   if (jte == jde) grid%xp%v1z(its:ite,jte,kts:kte) = 0.0

   !---------------------------------------------------------------------------
   ! [2.0] Perform forward FFT in x direction:
   !---------------------------------------------------------------------------

   ! [2.1] Apply (i',j',k -> i,j',k') transpose (v1z -> v1x).

   call da_transpose_z2x (grid)

   ! [2.2] Set up FFT parameters:
   
   idim = xbx%fft_ix
   jdim = grid%xp%jtex-grid%xp%jtsx+1

   vector_inc  = 1
   vector_jump = idim
   vector_size = idim - 1

   num_vectors = jdim

   work_area   = (vector_size+1)*num_vectors

   ! [2.3] Perform forward FFT:

   do k = grid%xp%ktsx, grid%xp%ktex
      ij = 0
      do j=grid%xp%jtsx, grid%xp%jtex
         do i=ids, ide
            ij=ij+1
            work_1d(ij) = grid%xp%v1x(i,j,k)
         end do

         do i=1, xbx%fft_pad_i
            ij=ij+1
            work_1d(ij) = 0.0
         end do
      end do

      call fft661(Forward_FFT, vector_inc, vector_jump, &
         num_vectors, vector_size, &
         xbx % fft_factors_x, xbx % trig_functs_x, &
         work_1d(1), work_area)
      ij = 0
      do j=grid%xp%jtsx, grid%xp%jtex
         do i=ids, ide
            ij=ij+1
            grid%xp%v1x(i,j,k) = work_1d(ij)
         end do

         do n=1, xbx%fft_pad_i
            i=(n-1)*xbx%pad_inc + 1
            ij=ij+1
            grid%xp%v2x(i,j,k) = work_1d(ij)
         end do
      end do
   end do

   !---------------------------------------------------------------------------
   ! [3.0] For each k-level, perform forward FFT in y direction, apply spectral 
   !        Poisson equation, and then perform inverse FFT in y direction:
   !---------------------------------------------------------------------------

   ! [3.1] Apply (i,j',k' -> i',j,k') transpose (v1x -> v1y).

   call da_transpose_x2y (grid)
   call da_transpose_x2y_v2 (grid)

   ! [3.2] Set up FFT parameters:

   idim = grid%xp%itey - grid%xp%itsy + 1 + xbx%pad_num
   jdim = xbx%fft_jy

   vector_inc  = idim
   vector_jump = 1
   vector_size = jdim - 1
   num_vectors = idim

   work_area   = (vector_size+1)*num_vectors

   ! [2.3] Perform forward FFT in j:

   do k = grid%xp%ktsy, grid%xp%ktey
      ij = 0
      do j=jds,jde
         do i=grid%xp%itsy, grid%xp%itey
            ij=ij+1
            work_1d(ij) = grid%xp%v1y(i,j,k)
         end do

         do n=1, xbx%pad_num
            i=xbx%pad_loc(n)
            ij=ij+1
            work_1d(ij) = grid%xp%v2y(i,j,k)
         end do
      end do

      do j=1, xbx%fft_pad_j
         do i=grid%xp%itsy, grid%xp%itey+xbx%pad_num
            ij=ij+1
            work_1d(ij) = 0.0
         end do
      end do

      call fft661(Forward_FFT, vector_inc, vector_jump, &
         num_vectors, vector_size, &
         xbx % fft_factors_y, xbx % trig_functs_y, &
         work_1d(1), work_area)

      !------------------------------------------------------------------------
      ! [4.0] Solve spectral Poisson equation:
      !------------------------------------------------------------------------

      ij = 0
      do j=jds, xbx%fft_jy
         do i=grid%xp%itsy, grid%xp%itey
            ij=ij+1
            work_1d(ij) = xbx%fft_coeffs(i,j)*work_1d(ij)
         end do
         do n=1, xbx%pad_num
            i=xbx%pad_pos(n)
            ij=ij+1
            work_1d(ij) = xbx%fft_coeffs(i,j)*work_1d(ij)
         end do
      end do

      ! [2.3] Reform gridpt. b using inverse double FST in i.

      call fft661(Inverse_FFT, vector_inc, vector_jump, &
         num_vectors, vector_size, &
         xbx % fft_factors_y, xbx % trig_functs_y, &
         work_1d(1), work_area)
                          
      ij = 0
      do j=jds, jde
         do i=grid%xp%itsy, grid%xp%itey
            ij=ij+1
            grid%xp%v1y(i,j,k) = work_1d(ij)
         end do

         do n=1, xbx%pad_num
            i=xbx%pad_loc(n)
            ij=ij+1
            grid%xp%v2y(i,j,k) = work_1d(ij)
         end do
      end do
   end do

   !---------------------------------------------------------------------------
   ! Perform inverse FFT in x direction:
   !---------------------------------------------------------------------------

   ! Apply (i',j,k' -> i,j',k') transpose (v1y -> v1x).

   call da_transpose_y2x (grid)
   call da_transpose_y2x_v2 (grid)

   ! Set up FFT parameters:
   
   idim = xbx%fft_ix
   jdim = grid%xp%jtex-grid%xp%jtsx+1

   vector_inc  = 1
   vector_jump = idim
   vector_size = idim - 1

   num_vectors = jdim

   work_area   = (vector_size+1)*num_vectors

   do k = grid%xp%ktsx, grid%xp%ktex
      ij = 0
      do j=grid%xp%jtsx, grid%xp%jtex
         do i=grid%xp%ids, grid%xp%ide
            ij=ij+1
            work_1d(ij) = grid%xp%v1x(i,j,k)
         end do

         do n=1, xbx%fft_pad_i
            i=(n-1)*xbx%pad_inc + 1
            ij=ij+1
            work_1d(ij) = grid%xp%v2x(i,j,k)
         end do
      end do

      call fft661(Inverse_FFT, vector_inc, vector_jump, &
         num_vectors, vector_size, &
         xbx % fft_factors_x, xbx % trig_functs_x, &
         work_1d(1), work_area)
      ij = 0
      do j=grid%xp%jtsx, grid%xp%jtex
         do i=ids, ide
            ij=ij+1
            grid%xp%v1x(i,j,k) = work_1d(ij)
         end do

         ij=ij+xbx%fft_pad_i
      end do
   end do

   ! Apply (i,j',k') -> i',j',k) transpose to restore v1z.

   call da_transpose_x2z (grid)

   !---------------------------------------------------------------------------
   ! [5.0] Tidy up:
   !---------------------------------------------------------------------------

   if (allocated(work_1d)) deallocate(work_1d)

   ! [2.5] Write data array into b:

   b(its:ite,jts:jte,kts:kte) = grid%xp%v1z(its:ite,jts:jte,kts:kte)

   if (trace_use) call da_trace_exit("da_solve_poissoneqn_fst")

end subroutine da_solve_poissoneqn_fst