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