da_vv_to_v_spectral.inc
References to this file elsewhere.
1 subroutine da_vv_to_v_spectral(ni, nj, max_wavenumber, inc, lenr, lensav, lenwrk, &
2 alp_size, r_cvsize, alp, wsave, int_wgts, rcv, field)
3
4 !-------------------------------------------------------------------------
5 ! Purpose: Performs gridpoint to spectral transformation on a sphere.
6 ! Note: Routine works for both regular and Gaussian latitude (latitude
7 ! integration weights contain necessary info).
8 !-------------------------------------------------------------------------
9
10 implicit none
11
12 integer, intent(in) :: ni ! Number of longitudes.
13 integer, intent(in) :: nj ! Number of latitudes.
14 integer, intent(in) :: r_cvsize ! Size of real control cv-array.
15 integer, intent(in) :: max_wavenumber ! Smallest scale required (ni/2 - 1).
16 integer, intent(in) :: inc ! Jump between elements of vector in array.
17 integer, intent(in) :: lenr ! FFT array dimension (at least inc*(n-1)+1).
18 real, intent(in) :: field(1:ni,1:nj) ! Gridpoint field.
19 real, intent(out) :: rcv(1:r_cvsize) ! Spectral modes.
20 integer, intent(in) :: lensav ! wsave dimension (n+int(log(real(ni)))+4).
21 integer, intent(in) :: lenwrk ! Dimension of work array.
22 integer, intent(in) :: alp_size ! Size of ALP vector.
23 real, intent(in) :: alp(1:alp_size) ! Associated Legendre Polynomials.
24 real, intent(in) :: wsave(1:lensav) ! Primes for FFT.
25 real, intent(in) :: int_wgts(1:nj) ! Legendre integration weights.
26
27 integer :: i, j, m, mm ! Loop counters.
28 integer :: sizec ! Size of complex cv-array
29 integer :: index_r, index_c ! Array index for complex v_fou
30 integer :: index_start, index_end ! Position markers in cv.
31 real :: r_fou(1:lenr) ! FFT array.
32 real :: work(1:lenwrk) ! FFT work array.
33 logical :: odd_longitudes
34 complex :: v_fou(1:nj,0:max_wavenumber)! Intermediate Fourier state.
35 complex :: r_leg(1:nj) ! Intermediate Fourier state.
36 complex, allocatable:: ccv(:) ! Spectral modes.
37
38 sizec = INT(0.5 * r_cvsize)
39 allocate (ccv(1:sizec))
40
41 if ((ni+1) / 2 == ni/2 + 1) then ! Odd number of longitudes:
42 odd_longitudes = .true.
43 else ! Even number of longitudes:
44 odd_longitudes = .false.
45 end if
46
47 !-------------------------------------------------------------------------
48 ! [1] Perform Adjoint of inverse Fourier decomposition in E-W direction:
49 !-------------------------------------------------------------------------
50
51 if ((ni+1) / 2 == ni/2 + 1) then ! Odd number of longitudes:
52 odd_longitudes = .true.
53 else ! Even number of longitudes:
54 odd_longitudes = .false.
55 end if
56
57 ! [1] Perform Fourier decomposition in E-W direction:
58
59 do j = 1, nj
60 r_fou(1:ni) = field(1:ni,j)
61 #ifdef FFTPACK
62 call rfft1f(ni, inc, r_fou, lenr, wsave, lensav, work, lenwrk, ierr)
63 #else
64 call da_error(__FILE__,__LINE__,(/"Needs to be compiled with FFTPACK"/))
65 #endif
66
67 if (odd_longitudes) then
68 v_fou(j,0) = CMPLX(r_fou(1), 0.0) ! m = 0 is real.
69 else
70 ! m = 0 is real, but pack R(NI/2) in imag m = 0:
71 v_fou(j,0) = CMPLX(r_fou(1), r_fou(ni))
72 end if
73
74 do m = 1, max_wavenumber
75 index_r = 2 * m
76 index_c = 2 * m + 1
77 v_fou(j,m) = CMPLX(r_fou(index_r), r_fou(index_c)) ! 2.0 * Fourier mode.
78 end do
79 end do
80
81 ! [2] Perform Legendre decomposition in N-S direction:
82
83 do m = 0, max_wavenumber
84 index_start = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1
85 index_end = index_start + max_wavenumber - m
86 r_leg(1:nj) = v_fou(1:nj,m)
87 call da_legtra (nj, max_wavenumber, alp_size, m, int_wgts, alp, r_leg, &
88 ccv(index_start:index_end))
89 end do
90
91 do i=1,sizec
92 mm = 2*i - 1
93 rcv(mm ) = real (ccv(i))
94 rcv(mm+1) = aimag(ccv(i))
95 end do
96 deallocate (ccv)
97
98 end subroutine da_vv_to_v_spectral
99
100