da_vtovv_spectral_adj.inc
References to this file elsewhere.
1 subroutine da_vtovv_spectral_adj(max_wavenumber, sizec, lenr, lenwrk, lensav, inc, &
2 alp_size, alp, wsave, power, rcv, field)
3
4 !----------------------------------------------------------------------
5 ! Purpose: Performs Adjoint of spectral to grid transformation on a sphere.
6 !----------------------------------------------------------------------
7
8 implicit none
9
10 integer, intent(in) :: max_wavenumber ! 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_wavenumber) ! Power spectrum
20 real, intent(out) :: rcv(1:2*sizec) ! Spectral modes.
21 real, intent(in) :: field(its:ite,jts:jte) ! Gridpoint field.
22
23 integer :: j, m, n ! Loop counters.
24 integer :: index_start ! Position markers in cv.
25 #ifdef DM_PARALLEL
26 integer :: index_end ! Position markers in cv.
27 #endif
28 integer :: index_r, index_c ! Array index for complex v_fou.
29
30 real :: r_fou(1:lenr) ! FFT array.
31 complex :: v_fou(its:ite,0:max_wavenumber)! Intermediate Fourier state.
32 complex :: ccv(1:sizec) ! Spectral modes.
33 complex :: ccv_local(1:sizec) ! Spectral modes.
34
35 integer :: l, js, je ! Loop counters.
36 integer :: index_m, index_j ! Markers.
37 complex :: sum_legtra ! Summation scalars.
38
39 integer :: jc, iequator, temp
40
41 #ifdef FFTPACK
42 real :: work(1:lenwrk) ! FFT work array.
43 #endif
44
45 if (trace_use) call da_trace_entry("da_vtovv_spectral_adj")
46
47 !----------------------------------------------------------------------
48 ! [1] Perform Adjoint of inverse Fourier decomposition in E-W direction:
49 !----------------------------------------------------------------------
50
51 v_fou = 0.0
52 do j = jts, jte
53 r_fou(its:ite) = field(its:ite,j)
54 #ifdef FFTPACK
55 call rfft1f(ide, inc, r_fou, lenr, wsave, lensav, work, lenwrk, ierr)
56 #else
57 call da_error(__FILE__,__LINE__,(/"Must compile with FFTPACK"/))
58 #endif
59
60 !----------------------------------------------------------------------
61 ! Adjust the output for adjoint test
62 !----------------------------------------------------------------------
63 r_fou = real(ite)/2.0 * r_fou
64 r_fou(its) = r_fou(its) * 2.0
65
66 ! if(.not. odd_longitudes) r_fou(ite) = 2.0*r_fou(ite)
67 ! make r_fou(ide) zero as there is no power computed for this wavenumber
68 r_fou(ite) = 0.0
69
70 v_fou(j,0) = CMPLX(r_fou(its), r_fou(ite))
71
72 do m = 1, max_wavenumber
73 index_r = 2 * m
74 index_c = 2 * m + 1
75 v_fou(j,m) = v_fou(j,m) + cmplx(r_fou(index_r),r_fou(index_c))
76 end do
77 end do
78
79 !----------------------------------------------------------------------
80 ! [2] Perform adjoint of inverse Legendre decomposition in N-S direction:
81 !----------------------------------------------------------------------
82
83 ccv_local(:) = 0.0
84
85 do m = 0, max_wavenumber
86 index_start = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 +1
87
88 index_m = m * (max_wavenumber + 1 - m) + m * (m + 1) / 2 + 1 - m
89
90 jc = (jde-jds+1)/2
91
92 iequator = mod(jde-jds+1, 2)
93
94 js = max(jts, jc+iequator+1)
95 je = min(jc+iequator, jte)
96
97 temp = (max_wavenumber + 1) * (max_wavenumber + 2) / 2
98
99 do l = m, max_wavenumber
100 sum_legtra = da_zero_complex
101
102 if (mod(l+m,2) == 1) then
103 do j = js, jte
104 index_j = (jds+jde - j - 1) * temp
105 sum_legtra = sum_legtra - v_fou(j,m) * alp(index_j + index_m + l)
106 end do
107 else
108 do j = js, jte
109 index_j = (jds+jde - j - 1) * temp
110 sum_legtra = sum_legtra + v_fou(j,m) * alp(index_j + index_m + l)
111 end do
112 end if
113
114 do j = jts, je
115 index_j = (j - 1) * temp
116 sum_legtra = sum_legtra + v_fou(j,m) * alp(index_j + index_m + l)
117 end do
118
119 ccv_local(index_start+l-m) = sum_legtra
120 end do
121 end do
122
123 #ifdef DM_PARALLEL
124 index_start = 1
125 index_end = max_wavenumber + &
126 max_wavenumber * (max_wavenumber + 1) / 2 + 1
127
128 n = index_end - index_start + 1
129 call mpi_allreduce(ccv_local(index_start:index_end), &
130 ccv(index_start:index_end), n, true_mpi_complex, mpi_sum, &
131 comm, ierr)
132 #else
133 ccv(:) = ccv_local(:)
134 #endif
135
136 !----------------------------------------------------------------------
137 ! [2] Apply Power spectrum
138 !-------------------------------------------------------------------------
139
140 if (.not. test_transforms) call da_apply_power(power, max_wavenumber, ccv, sizec)
141
142 do n = 1, sizec
143 rcv(2*n - 1) = real (ccv(n))
144 rcv(2*n ) = aimag(ccv(n))
145 end do
146
147 if (trace_use) call da_trace_exit("da_vtovv_spectral_adj")
148
149 end subroutine da_vtovv_spectral_adj
150
151