da_vtovv_spectral.inc
References to this file elsewhere.
1 subroutine da_vtovv_spectral(max_wave, sizec, lenr, lenwrk, lensav, inc, &
2 alp_size, alp, wsave, power, rcv, field)
3
4 !-------------------------------------------------------------------------
5 ! Purpose: Performs spectral to gridpoint transformation on a sphere.
6 !----------------------------------------------------------------------
7
8 implicit none
9
10 integer, intent(in) :: max_wave ! Max total wavenumber.
11 integer, intent(in) :: sizec ! Size of packed spectral array.
12 integer, intent(in) :: lenr ! FFT info.
13 integer, intent(in) :: lenwrk ! FFT info.
14 integer, intent(in) :: lensav ! FFT info.
15 integer, intent(in) :: inc ! FFT info.
16 integer, intent(in) :: alp_size ! Size of alp array.
17 real, intent(in) :: alp(1:alp_size) ! Associated Legendre Polynomials.
18 real, intent(in) :: wsave(1:lensav) ! Primes for FFT.
19 real, intent(in) :: power(0:max_wave) ! Power spectrum
20 real, intent(in) :: rcv(1:2*sizec) ! Spectral modes.
21 real, intent(out) :: field(its:ite,jts:jte) ! Gridpoint field.
22
23 integer :: j, l,m, n ! Loop counters.
24 integer :: index_start ! Position markers in cv.
25 integer :: index_r, index_c ! Array index for complex v_fou.
26
27 real :: r_fou(1:lenr) ! FFT array.
28 complex :: v_fou(its:ite,0:max_wave)! Intermediate Fourier state.
29 complex :: ccv(1:sizec) ! Spectral modes.
30 integer :: index_m, index_j
31 integer :: jc, js, je, iequator
32 complex :: sum_legtra ! Summation scalars.
33
34 #ifdef FFTPACK
35 real :: work(1:lenwrk) ! FFT work array.
36 #endif
37
38
39 if (trace_use) call da_trace_entry("da_vtovv_spectral")
40
41 !----------------------------------------------------------------------------
42 ! [1] Create complex array from read array:
43 !----------------------------------------------------------------------------
44
45 v_fou = 0.0
46 do n = 1, sizec
47 ccv(n) = CMPLX(rcv(2*n-1), rcv(2*n))
48 end do
49
50 !----------------------------------------------------------------------------
51 ! [2] Apply power spectrum
52 !----------------------------------------------------------------------------
53
54 if (.not. test_transforms) call da_apply_power(power, max_wave, ccv, sizec)
55
56 !----------------------------------------------------------------------------
57 ! [3] Perform inverse Legendre decomposition in N-S direction:
58 !----------------------------------------------------------------------------
59
60 do m = 0, max_wave
61 index_start = m * (max_wave + 1 - m) + m * (m + 1) / 2 + 1
62
63 index_m = m * (max_wave + 1 - m) + m * (m + 1) / 2 + 1 - m
64
65 jc = (jde-jds+1)/2
66
67 iequator = mod(jde-jds+1, 2)
68
69 je = min(jc+iequator, jte)
70
71 do j = jts, je
72 index_j = (j - 1) * (max_wave + 1) * (max_wave + 2) / 2
73
74 v_fou(j,m) = sum(ccv(index_start:index_start-m+max_wave) * &
75 alp(index_j+index_m+m:index_j+index_m+max_wave))
76 end do
77
78 js = max(jts, jc+iequator+1)
79
80 do j = js, jte
81 index_j = (jds+jde - j - 1) * (max_wave + 1) * (max_wave + 2) / 2
82
83 sum_legtra = da_zero_complex
84 do l = m, max_wave
85 ! Calculate second quadrant values:
86 if(mod(l+m,2) == 1) then
87 sum_legtra = sum_legtra - ccv(index_start-m+l) * alp(index_j + index_m + l)
88 else
89 sum_legtra = sum_legtra + ccv(index_start-m+l) * alp(index_j + index_m + l)
90 end if
91 end do
92 v_fou(j,m) = sum_legtra
93 end do
94 end do
95
96 !----------------------------------------------------------------------------
97 ! [4] Perform inverse Fourier decomposition in E-W direction:
98 !----------------------------------------------------------------------------
99
100 do j = jts, jte
101 r_fou(its) = real(v_fou(j,0)) ! R(m=0) is real.
102 ! r_fou(ite) = aimag(v_fou(j,0)) ! R(m=NI/2) is real, but packed in imag m = 0)
103 ! make r_fou(ide) zero as there is no power computed corresponding to this wavenumber
104 r_fou(ite) = 0.0
105
106 do m = 1, max_wave
107 index_r = 2 * m
108 index_c = 2 * m + 1
109 r_fou(index_r) = real(v_fou(j,m))
110 r_fou(index_c) = aimag(v_fou(j,m))
111 end do
112
113 #ifdef FFTPACK
114 call rfft1b(ide, inc, r_fou, lenr, wsave, lensav, work, lenwrk, ierr)
115 #else
116 call da_error(__FILE__,__LINE__,(/"Must compile with FFTPACK"/))
117 #endif
118 field(its:ite,j) = r_fou(its:ite)
119 end do
120
121 if (trace_use) call da_trace_exit("da_vtovv_spectral")
122
123 end subroutine da_vtovv_spectral
124
125