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