subroutine da_calc_power_spectrum(max_wave, power)

   !-----------------------------------------------------------------------
   ! Purpose: Calculate power spectrum
   !-----------------------------------------------------------------------
    
   implicit none
    
   integer, parameter  :: nj = 200          ! #Gaussian lats (even).

   integer, intent(in) :: max_wave          ! Smallest wavenumber for grid.
   real*8, intent(out) :: power(0:max_wave) ! Power spectrum.

   real                :: glats(1:nj)       ! Gaussian latitudes.
   real                :: gwgts(1:nj)       ! Gaussian weights.
   real                :: sinlat(1:nj)      ! sine.
   real                :: coslat(1:nj)      ! cosine.

   integer             :: max_wave_nj       ! Maximum wavenumber.
   integer             :: i,j, l            ! Loop counters.
   real                :: corr_scale        ! Correlation scale.
   real                :: corr_scale_inv    ! 1/corr_scale
   real                :: variance          ! Variance = sum(power).
   real                :: d(1:nj)           ! Temp array.
   real                :: corr(1:nj)        ! Correlation function.
   logical             :: odd               ! true if odd.
   real, allocatable   :: alp(:,:)          ! Associated Legendre Polynomials.
   real, allocatable   :: power_nj(:)       ! Power spectrum.
   character (len=filename_len) :: filename
   
   do i=1,num_alpha_corr_types
      call da_get_unit(alpha_corr_unit1(i))
      call da_get_unit(alpha_corr_unit2(i))
      write (unit=filename,fmt='(A,I1)') "alpha_corr1_",i
      open(unit=alpha_corr_unit1(i),file=filename,status="replace")
      write (unit=filename,fmt='(A,I1)') "alpha_corr2_",i
      open(unit=alpha_corr_unit2(i),file=filename,status="replace")
   end do

   !-----------------------------------------------------------------------------
   ! [1] Switch lats from -pi/2 to pi/2 to 0 to pi:
   !-----------------------------------------------------------------------------

   max_wave_nj = nj / 2 - 1
   allocate(alp(1:nj,0:max_wave_nj))
   allocate(power_nj(0:max_wave_nj))

   call da_get_gausslats(nj, glats, gwgts, sinlat, coslat)
   glats = glats + 0.5 * pi

   ! Get m=0 Associated Legendre Polynomials:

   do l = 0, max_wave_nj
      odd = .false.
      if (real(int(0.5 * real(l))) /= 0.5 * l) odd = .true.

      do j = 1, nj
         call da_asslegpol(l, 0, sinlat(j), coslat(j), alp(j,l))
         ! Reverse order of alps to account for latitude/angle difference:
         if (odd) alp(j,l) = -alp(j,l)
      end do
   end do

   !-----------------------------------------------------------------------------
   ! [2] Define correlation function:
   !-----------------------------------------------------------------------------

   corr_scale = alpha_corr_scale / earth_radius 
   corr_scale_inv = 1.0 / corr_scale

   do j = 1, nj
      ! d(j) = 0.5 * glats(j) * corr_scale_inv
      d(j) = glats(j) * corr_scale_inv

      if (alpha_corr_type == alpha_corr_type_exp) then ! Exponential.
         corr(j) = exp(-d(j))
      else if (alpha_corr_type == alpha_corr_type_soar) then ! SOAR
         d(j) = 2.0 * d(j)
         corr(j) = (1.0 + d(j)) * exp(-d(j))
      else if (alpha_corr_type == alpha_corr_type_gaussian) then ! Gaussian
         corr(j) = exp(-d(j) * d(j))
      end if
   end do

   do j = 1, nj
      write(unit=alpha_corr_unit1(alpha_corr_type),fmt='(i4,2f12.4)') &
        j, earth_radius * glats(j), corr(j)
   end do

   !--------------------------------------------------------------------------
   ! [3] Calculate power spectra:
   !--------------------------------------------------------------------------

   ! Calculate power spectrum (and truncate if has -ve values).
   ! Power spectrum at this stage is is the Dl=sqrt(2l+1)*Bl of Boer(1983).
   ! This ensures the total variance = sum(Dl).

   power_nj(:) = 0.0
   do l = 0, max_wave_nj
      power_nj(l) = 0.5 * sqrt(2.0 * real(l) + 1.0) * &
                    sum(gwgts(1:nj) * corr(1:nj) * alp(1:nj,l))

      if (power_nj(l) < 0.0) power_nj(l) = 0.0
   end do
   write(unit=stdout,fmt='(a,2f12.5)')' Total, unscale variance = ', &
                sum(power_nj(0:max_wave_nj))

   ! Rescale so variance = 1 (take out later?):
   variance = sum(power_nj(0:max_wave_nj))
   power_nj(0:max_wave_nj) = power_nj(0:max_wave_nj) / variance

   do l = 0, max_wave_nj
      write(unit=alpha_corr_unit2(alpha_corr_type),fmt='(i4,2f12.4)') &
        l, power_nj(l), sum(power_nj(0:l))
   end do

   write(unit=stdout,fmt='(a,2i6)')' Total, truncated wave_max = ', &
                     max_wave_nj, max_wave
   write(unit=stdout,fmt='(a,2f12.5)')' Total, truncated variance = ', &
                sum(power_nj(0:max_wave_nj)), sum(power_nj(0:max_wave))

   power(0:max_wave) = power_nj(0:max_wave)

   ! Add compactly supported correlation from calc_globalspectral later?

   deallocate(alp)
   deallocate(power_nj)
   
   do i=1,num_alpha_corr_types
      close (alpha_corr_unit1(i))
      close (alpha_corr_unit2(i))
      call da_free_unit (alpha_corr_unit1(i))
      call da_free_unit (alpha_corr_unit2(i))
   end do

end subroutine da_calc_power_spectrum