subroutine da_vtovv_spectral_adj(max_wavenumber, sizec, lenr, lenwrk, lensav, inc, &
                                    alp_size, alp, wsave, power, rcv, field)

   !----------------------------------------------------------------------
   ! Purpose: Performs Adjoint of spectral to grid transformation on a sphere.
   !----------------------------------------------------------------------

   implicit none

   integer, intent(in) :: max_wavenumber             ! Max total wavenumber.
   integer, intent(in) :: sizec                  ! Size of packed spectral array.
   integer, intent(in) :: lenr                       ! FFT info.
   integer, intent(in) :: lenwrk                     ! FFT info.
   integer, intent(in) :: lensav                     ! FFT info.
   integer, intent(in) :: inc                        ! FFT info.
   integer, intent(in) :: alp_size                   ! Size of alp array.
   real, intent(in)    :: alp(1:alp_size)            ! Associated Legendre Polynomials.
   real, intent(in)    :: wsave(1:lensav)            ! Primes for FFT.
   real*8, intent(in)    :: power(0:max_wavenumber)    ! Power spectrum
   real, intent(out)   :: rcv(1:2*sizec)             ! Spectral modes.
   real, intent(in)    :: field(its:ite,jts:jte)     ! Gridpoint field.

   integer             :: j, m, n                    ! Loop counters.
   integer             :: index_start                ! Position markers in cv.
#ifdef DM_PARALLEL
   integer             :: index_end                  ! Position markers in cv.
#endif
   integer             :: index_r, index_c           ! Array index for complex v_fou.

   real                :: r_fou(1:lenr)              ! FFT array.
   complex             :: v_fou(its:ite,0:max_wavenumber)! Intermediate Fourier state.
   complex             :: ccv(1:sizec)               ! Spectral modes.
   complex             :: ccv_local(1:sizec)         ! Spectral modes.

   integer              :: l, js, je           ! Loop counters.
   integer              :: index_m, index_j    ! Markers.
   complex              :: sum_legtra          ! Summation scalars.

   integer              :: jc, iequator, temp

#ifdef FFTPACK
   real                :: work(1:lenwrk)             ! FFT work array.
#endif

   if (trace_use) call da_trace_entry("da_vtovv_spectral_adj")

   !----------------------------------------------------------------------
   ! [1] Perform Adjoint of inverse Fourier decomposition in E-W direction:
   !----------------------------------------------------------------------

   v_fou = 0.0
   do j = jts, jte
      r_fou(its:ite) = field(its:ite,j) 
#ifdef FFTPACK
      call rfft1f(ide, inc, r_fou, lenr, wsave, lensav, work, lenwrk, ierr)
#else
      call da_error(__FILE__,__LINE__,(/"Must compile with FFTPACK"/))
#endif

      !----------------------------------------------------------------------
      ! Adjust the output for adjoint test
      !----------------------------------------------------------------------
      r_fou      =  real(ite)/2.0 * r_fou
      r_fou(its) =  r_fou(its)   * 2.0       

      ! if(.not. odd_longitudes) r_fou(ite) = 2.0*r_fou(ite)   
      ! make r_fou(ide) zero as there is no power computed for this wavenumber
      r_fou(ite) = 0.0

      v_fou(j,0) = CMPLX(r_fou(its), r_fou(ite))

      do m = 1, max_wavenumber
         index_r = 2 * m
         index_c = 2 * m + 1
         v_fou(j,m)  = v_fou(j,m) + cmplx(r_fou(index_r),r_fou(index_c))
      end do
   end do

   !----------------------------------------------------------------------
   ! [2] Perform adjoint of inverse Legendre decomposition in N-S direction:
   !----------------------------------------------------------------------

   ccv_local(:) = 0.0

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

      index_m     = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 - m

      jc = (jde-jds+1)/2

      iequator = mod(jde-jds+1, 2)

      js = max(jts, jc+iequator+1)
      je = min(jc+iequator, jte)

      temp = (max_wavenumber + 1) * (max_wavenumber + 2) / 2

      do l = m, max_wavenumber
         sum_legtra = da_zero_complex

         if (mod(l+m,2) == 1) then
            do j = js, jte
               index_j = (jds+jde - j - 1) * temp
               sum_legtra = sum_legtra - v_fou(j,m) * alp(index_j + index_m + l)
            end do
         else
            do j = js, jte
               index_j = (jds+jde - j - 1) * temp
               sum_legtra = sum_legtra + v_fou(j,m) * alp(index_j + index_m + l)
            end do
         end if

         do j = jts, je
            index_j = (j - 1) * temp
            sum_legtra = sum_legtra + v_fou(j,m) * alp(index_j + index_m + l) 
         end do
   
         ccv_local(index_start+l-m) = sum_legtra
      end do
   end do

#ifdef DM_PARALLEL
   index_start = 1
   index_end   = max_wavenumber + &
      max_wavenumber * (max_wavenumber + 1) / 2 + 1

   n = index_end - index_start + 1
   call mpi_allreduce(ccv_local(index_start:index_end), &
                      ccv(index_start:index_end), n, true_mpi_complex, mpi_sum, &
                      comm, ierr)
#else
   ccv(:) = ccv_local(:)                      
#endif

   !----------------------------------------------------------------------
   ! [2] Apply Power spectrum
   !-------------------------------------------------------------------------

   if (.not. test_transforms) call da_apply_power(power, max_wavenumber, ccv, sizec)

   do n = 1, sizec
      rcv(2*n - 1) = real (ccv(n))
      rcv(2*n    ) = aimag(ccv(n))
   end do 

   if (trace_use) call da_trace_exit("da_vtovv_spectral_adj")

end subroutine da_vtovv_spectral_adj