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