da_solve_poissoneqn_fst_adj.inc
References to this file elsewhere.
1 subroutine da_solve_poissoneqn_fst_adj(xbx, xp, b, del2b)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Solve Del**2 B = A for B with zero field boundary conditions.
5 !
6 ! Method: 1) Compute spectral del2b using double forward FST.
7 ! 2) Calculate spectral b.
8 ! 3) Reform gridpt. b using inverse double FST.
9 ! Note no mean removed (would make b.ne.0 at boundaries) so
10 ! arbitrary constant still present.
11 !---------------------------------------------------------------------------
12
13 implicit none
14
15 type (xbx_type),intent(in) :: xbx ! Header & non-gridded vars.
16 type (xpose_type),intent(inout) :: xp ! Dimensions and xpose buffers.
17 real, intent(in) :: b(ims:ime,jms:jme,kms:kme) ! B.
18 real, intent(out) :: del2b(ims:ime,jms:jme,kms:kme) ! Del**2 B.
19
20 integer :: vector_inc ! Increment between FST data.
21 integer :: vector_jump ! Jump between start of vectors.
22 integer :: num_vectors ! Number of FSTs to perform.
23 integer :: vector_size ! Of form 2**p 3**q 5**r for FSTs.
24 integer :: work_area ! Dimension for FST routine.
25 integer :: idim ! Size of 1st dimension for FST.
26 integer :: jdim ! Size of 2nd dimension for FST.
27
28 integer :: i, j, k, n, ij ! loop counter
29
30 real, allocatable :: work_1d(:) ! FFT work array for x dir.
31
32 !---------------------------------------------------------------------------
33 ! [1.0] Initialise:
34 !---------------------------------------------------------------------------
35
36 ! Calculate work space needed.
37
38 n = max(xbx%fft_ix*(xp%jtex-xp%jtsx+1), &
39 xbx%fft_jy*(xp%itey-xp%itsy+1+xbx%pad_num))
40
41 ! Allocate work arrays.
42 allocate(work_1d(1:n))
43
44 ! Copy b for transpose.
45
46 xp%v1z(its:ite,jts:jte,kts:kte) = b(its:ite,jts:jte,kts:kte)
47
48 !---------------------------------------------------------------------------
49 ! [2.0] Perform forward FFT in x direction:
50 !---------------------------------------------------------------------------
51
52 ! [2.1] Apply (i',j',k -> i,j',k') transpose (v1z -> v1x).
53
54 call da_transpose_z2x (xp)
55
56 ! [2.2] Set up FFT parameters:
57
58 idim = xbx%fft_ix
59 jdim = xp%jtex - xp%jtsx + 1
60
61 vector_inc = 1
62 vector_jump = idim
63 vector_size = idim - 1
64
65 num_vectors = jdim
66
67 work_area = (vector_size+1)*num_vectors
68
69 ! [2.3] Perform forward FFT:
70
71 do k = xp%ktsx, xp%ktex
72 ij = 0
73 do j=xp%jtsx, xp%jtex
74 do i=xp%ids, xp%ide
75 ij=ij+1
76 work_1d(ij) = xp%v1x(i,j,k)
77 end do
78
79 do n=1, xbx%fft_pad_i
80 ij=ij+1
81 work_1d(ij) = 0.0
82 end do
83 end do
84
85 call fft661(Forward_FFT, vector_inc, vector_jump, &
86 num_vectors, vector_size, &
87 xbx % fft_factors_x, xbx % trig_functs_x, &
88 work_1d(1), work_area)
89 ij = 0
90 do j=xp%jtsx, xp%jtex
91 do i=xp%ids, xp%ide
92 ij=ij+1
93 xp%v1x(i,j,k) = work_1d(ij)
94 end do
95
96 do n=1, xbx%fft_pad_i
97 i=(n-1)*xbx%pad_inc+1
98 ij=ij+1
99 xp%v2x(i,j,k) = work_1d(ij)
100 end do
101 end do
102 end do
103
104 !---------------------------------------------------------------------------
105 ! [3.0] For each k-level, perform forward FFT in y direction, apply spectral
106 ! Poisson equation, and then perform inverse FFT in x direction:
107 !---------------------------------------------------------------------------
108
109 ! [3.1] Apply (i,j',k' -> i',j,k') transpose (v1x -> v1y).
110
111 call da_transpose_x2y (xp)
112 call da_transpose_x2y_v2 (xp)
113
114 ! [3.2] Set up FFT parameters:
115
116 idim = xp%itey - xp%itsy + 1 + xbx%pad_num
117 jdim = xbx%fft_jy
118
119 vector_inc = idim
120 vector_jump = 1
121 vector_size = jdim - 1
122 num_vectors = idim
123
124 work_area = (vector_size+1)*num_vectors
125
126 do k = xp%ktsy, xp%ktey
127 ij = 0
128 do j=xp%jds, xp%jde
129 do i=xp%itsy, xp%itey
130 ij=ij+1
131 work_1d(ij) = xp%v1y(i,j,k)
132 end do
133
134 do n=1, xbx%pad_num
135 i=xbx%pad_loc(n)
136 ij=ij+1
137 work_1d(ij) = xp%v2y(i,j,k)
138 end do
139 end do
140
141 do j=1, xbx%fft_pad_j
142 do i=xp%itsy, xp%itey+xbx%pad_num
143 ij=ij+1
144 work_1d(ij) = 0.0
145 end do
146 end do
147
148 call fft661(Forward_FFT, vector_inc, vector_jump, &
149 num_vectors, vector_size, &
150 xbx % fft_factors_y, xbx % trig_functs_y, &
151 work_1d(1), work_area)
152
153 !------------------------------------------------------------------------
154 ! [4.0] Solve spectral Poisson equation:
155 !------------------------------------------------------------------------
156
157 ij = 0
158 do j=xp%jds, xbx%fft_jy
159 do i=xp%itsy, xp%itey
160 ij=ij+1
161 work_1d(ij) = xbx%fft_coeffs(i,j)*work_1d(ij)/xbx%fft_adjoint_factor
162 end do
163 do n=1, xbx%pad_num
164 i=xbx%pad_pos(n)
165 ij=ij+1
166 work_1d(ij) = xbx%fft_coeffs(i,j)*work_1d(ij)/xbx%fft_adjoint_factor
167 end do
168 end do
169
170 ! Reform gridpt. b using inverse double FST in i.
171
172 call fft661(Inverse_FFT, vector_inc, vector_jump, &
173 num_vectors, vector_size, &
174 xbx % fft_factors_y, xbx % trig_functs_y, &
175 work_1d(1), work_area)
176 ij = 0
177 do j=xp%jds, xp%jde
178 do i=xp%itsy, xp%itey
179 ij=ij+1
180 xp%v1y(i,j,k) = work_1d(ij)
181 end do
182
183 do n=1, xbx%pad_num
184 i=xbx%pad_loc(n)
185 ij=ij+1
186 xp%v2y(i,j,k) = work_1d(ij)
187 end do
188 end do
189 end do
190
191 !---------------------------------------------------------------------------
192 ! Perform inverse FFT in x direction:
193 !---------------------------------------------------------------------------
194
195 ! Apply (i',j,k' -> i,j',k') transpose (v1y -> v1x).
196
197 call da_transpose_y2x (xp)
198 call da_transpose_y2x_v2 (xp)
199
200 ! Set up FFT parameters:
201
202 idim = xbx%fft_ix
203 jdim = xp%jtex - xp%jtsx + 1
204
205 vector_inc = 1
206 vector_jump = idim
207 vector_size = idim - 1
208
209 num_vectors = jdim
210
211 work_area = (vector_size+1)*num_vectors
212
213 ! [2.3] Perform forward FFT:
214
215 do k = xp%ktsx, xp%ktex
216 ij = 0
217 do j=xp%jtsx, xp%jtex
218 do i=xp%ids, xp%ide
219 ij=ij+1
220 work_1d(ij) = xp%v1x(i,j,k)
221 end do
222
223 do n=1, xbx%fft_pad_i
224 i=(n-1)*xbx%pad_inc+1
225 ij=ij+1
226 work_1d(ij) = xp%v2x(i,j,k)
227 end do
228 end do
229
230 call fft661(Inverse_FFT, vector_inc, vector_jump, &
231 num_vectors, vector_size, &
232 xbx % fft_factors_x, xbx % trig_functs_x, &
233 work_1d(1), work_area)
234 ij = 0
235 do j=xp%jtsx, xp%jtex
236 do i=xp%ids, xp%ide
237 ij=ij+1
238 xp%v1x(i,j,k) = work_1d(ij)*xbx%fft_adjoint_factor
239 end do
240
241 ij=ij+xbx%fft_pad_i
242 end do
243 end do
244
245 ! Apply (i,j',k') -> i',j',k) transpose to restore v1z.
246
247 call da_transpose_x2z (xp)
248
249 call wrf_dm_halo(xp%domdesc,xp%comms,xp%halo_id2)
250
251 ! [2.5] Write data array into del2b:
252
253 del2b(ims:ime,jms:jme,kms:kme) = xp%v1z(ims:ime,jms:jme,kms:kme)
254
255 !---------------------------------------------------------------------------
256 ! [5.0] Tidy up:
257 !---------------------------------------------------------------------------
258
259 if (allocated(work_1d)) deallocate(work_1d)
260
261 end subroutine da_solve_poissoneqn_fst_adj
262
263