subroutine da_test_spectral (max_wave, sizec, xbx, field)

   !-----------------------------------------------------------------------
   ! Purpose: TBD
   !-----------------------------------------------------------------------

   implicit none

   integer,         intent(in) :: max_wave ! Maximum wavenumber.
   integer,         intent(in) :: sizec    ! Size of complex cv.
   type (xbx_type), intent(in) :: xbx  ! For header & non-grid arrays.
   real,            intent(in) :: field(its:ite,jts:jte)    ! Gridpoint field.

   real    :: field_out(its:ite,jts:jte)
   real*8  :: power(0:max_wave)   ! Power spectrum
   real    :: rcv(1:2*sizec)      ! Spectral modes.
   real    :: rcv_out(1:2*sizec)  ! Spectral modes.
   integer :: m,mm, index_start, index_end
   complex :: r_leg(1:jde)     
   complex :: ccv(1:sizec)     ! Spectral modes.
   complex :: ccv1(1:sizec)    ! Spectral modes.
   real    :: den, num, xx

   if (trace_use) call da_trace_entry("da_test_spectral")

   write (unit=stdout, fmt='(A)') &
      ' Test orthogonality of Associated Legendre Polynomials:'

   ! Initialise Power spectrum
   power  = 1.0
   call da_setlegpol_test( jde, max_wave, xbx%alp_size, xbx%int_wgts, xbx%alp )

   write(unit=stdout,fmt='(A)') &
      ' Test invertibility of spectral (Fourier, Legendre) transforms:'

   ! Gridpoint to spectral:
   rcv = 0.0
   call da_vv_to_v_spectral( xbx%ni, xbx%nj, max_wave, xbx%inc, xbx%lenr, xbx%lensav, xbx%lenwrk, &
                             xbx%alp_size, 2*sizec, xbx%alp, xbx%wsave, xbx%int_wgts, rcv, field)

   field_out = 0.0
   ! Spectral to gridpoint:
   call da_vtovv_spectral( max_wave, sizec, &
                             xbx % lenr, xbx % lenwrk, xbx % lensav, &
                             xbx % inc, xbx % alp_size, xbx % alp, &
                             xbx % wsave, power, rcv, field_out)

   write(unit=stdout,fmt='(1x,a,e30.10)') &
      'Domain-Averaged (Grid->Spectral->Grid) Error = ', &
              sqrt( sum( ( field_out(1:xbx%ni,1:xbx%nj) - field(1:xbx%ni,1:xbx%nj) )**2 ) / &
                    sum( field(1:xbx%ni,1:xbx%nj)**2 ) )
   rcv_out = 0.0
   
   ! Gridpoint to spectral (again):
   call da_vv_to_v_spectral( xbx%ni, xbx%nj, max_wave, xbx%inc, xbx%lenr, xbx%lensav, xbx%lenwrk, &
                      xbx%alp_size, 2*sizec, xbx%alp, xbx%wsave, xbx%int_wgts, rcv_out, field_out)

   rcv_out(1:2*sizec) = rcv_out(1:2*sizec) - rcv(1:2*sizec) ! Create difference for test diags.
    
   write(unit=stdout,fmt='(1x,a,e30.10)') &
      ' Domain-Averaged (Spectral->Grid->Spectral) Error = ', &
                       sqrt( sum( rcv_out(1:2*sizec)**2 ) ) / sqrt( sum( rcv(1:2*sizec)**2 ) )

   ! Adjoint test for Spectral Transform
   rcv_out = 0.0
   call da_vtovv_spectral_adj( max_wave, sizec, &
                                 xbx % lenr, xbx % lenwrk, xbx % lensav, &
                                 xbx % inc, xbx % alp_size, xbx % alp, &
                                 xbx % wsave, power, rcv_out, field_out)

   write(unit=stdout,fmt='(A)') ' Adjoint test result for  Spectral -> Grid : '
   write(unit=stdout,fmt='(1x,a,e30.10)') &
      ' LHS  ( LX.LX)       = ',&
      sum( field_out(1:xbx%ni,1:xbx%nj)*field_out(1:xbx%ni,1:xbx%nj) ) 
   write(unit=stdout,fmt='(1x,a,e30.10)') &
      ' RHS  (  X.L^TLX )   = ', sum( rcv(1:2*sizec)*rcv_out(1:2*sizec) ) 

   ! Inverse test for Legendre Transform

   write(unit=stdout,fmt='(A)') '  Inverse and Adjoint Legendre test result:'

   do m = 0, max_wave
      index_start = m * ( max_wave + 1 - m ) + m * ( m + 1 ) / 2 + 1
      index_end   = index_start + max_wave - m

      do mm = index_start, index_end
         if (2*mm > 2*sizec) then
            call da_error(__FILE__,__LINE__, &
               (/"rcv_out index bounce"/))
         end if
         ccv(mm) = cmplx (rcv_out(2*mm-1), rcv_out(2*mm))
      end do
      r_leg = 0.0
      call da_legtra_inv( jds, jde, jts, jte, max_wave, xbx%alp_size, m, xbx%alp, &
                          ccv(index_start:index_end), r_leg )

      ccv1(index_start:index_end) = 0.0
      call da_legtra ( xbx%nj, max_wave, xbx%alp_size, m, xbx%int_wgts, xbx%alp, r_leg, &
                       ccv1(index_start:index_end) )
      ccv1(index_start:index_end) = ccv1(index_start:index_end) - &
                                    ccv(index_start:index_end) 
      num =  sum ( real(ccv1(index_start:index_end))*real(ccv1(index_start:index_end))+&
         aimag(ccv1(index_start:index_end))* aimag(ccv1(index_start:index_end)) )     
      den =  sum ( real(ccv(index_start:index_end))*real(ccv(index_start:index_end))+&
         aimag(ccv(index_start:index_end))* aimag(ccv(index_start:index_end)) )     
      write(unit=stdout,fmt='(A,I4,A,E30.20)') &
         'For zonal wave number',m,' difference ',sqrt(num)/sqrt(den)

      xx = sum( real(r_leg(1:xbx%nj))* real(r_leg(1:xbx%nj))+ &
               aimag(r_leg(1:xbx%nj))*aimag(r_leg(1:xbx%nj)) )
      write(unit=stdout,fmt='(a,i5,a,e30.20)') 'For Wave = ',m,' LX.LX    = ',xx

      ccv1(index_start:index_end) = 0.0
      call da_legtra_inv_adj( jds, jde, jts, jte, max_wave, xbx%alp_size, m, xbx%alp, &
                              ccv1(index_start:index_end), r_leg )
      xx = sum( real(ccv(index_start:index_end))*     &
                real(ccv1(index_start:index_end))     +&
               aimag(ccv(index_start:index_end))* &
               aimag(ccv1(index_start:index_end)) )
      write(unit=stdout,fmt='(a,i5,a,e30.20,/)') 'For Wave = ',m,' X.L^TLX  = ',xx
   end do

   if (trace_use) call da_trace_exit("da_test_spectral")

end subroutine da_test_spectral