var/da/da_spectral/da_calc_power.inc

   1 subroutine da_calc_power(max_wavenumber, r_cvsize, rcv, power)
   2
   3    !----------------------------------------------------------------------
   4    ! Purpose: Performs spectral to gridpoint transformation on a sphere.
   5    !----------------------------------------------------------------------
   6
   7    implicit none
   8
   9    integer, intent(in) :: max_wavenumber ! Smallest scale required (ni/2 - 1)
  10    integer, intent(in) :: r_cvsize               ! Size of rcv vector.
  11    real, intent(in)    :: rcv(1:r_cvsize)        ! Spectral modes.
  12    real, intent(out)   :: power(0:max_wavenumber)! Power spectrum.
  13
  14    integer             :: n, m                   ! Loop counters.
  15    integer             :: index, indexx          ! Position markers in cv.
  16    integer             :: cv_size                ! Complex cv size.
  17    complex, allocatable:: ccv(:)                 ! Complex control variable.
  18
  19    ! Create complex array:
  20    cv_size = r_cvsize / 2
  21    allocate(ccv(1:cv_size))
  22    do index = 1, cv_size
  23       indexx = 2 * index - 1
  24       ccv(index) = cmplx(rcv(indexx), rcv(indexx+1))
  25    end do
  26
  27    power(:) = 0.0
  28
  29    ! Calculate power spectrum from input 1D spectral mode vector:
  30
  31    do n = 0, max_wavenumber
  32
  33       ! First consider m=0:
  34       m = 0
  35       index = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 + n - m
  36       power(n) = real(ccv(index))**2
  37
  38       ! Now add m>0 terms:
  39       do m = 1, n
  40          index = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 + n - m
  41          power(n) = power(n) + 2.0 *(real(ccv(index))**2 + aimag(ccv(index))**2)
  42       end do
  43
  44    end do
  45
  46    deallocate(ccv)
  47
  48 end subroutine da_calc_power
  49
  50