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