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