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