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