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