Home

Context Navigation

source: palm/trunk/SOURCE/poisfft.f90 @ 1385

Last change on this file since 1385 was 1321, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 44.8 KB

Line
1	MODULE poisfft_mod
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2014 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: poisfft.f90 1321 2014-03-20 09:40:40Z raasch $
27	!
28	! 1320 2014-03-20 08:40:49Z raasch
29	! ONLY-attribute added to USE-statements,
30	! kind-parameters added to all INTEGER and REAL declaration statements,
31	! kinds are defined in new module kinds,
32	! old module precision_kind is removed,
33	! revision history before 2012 removed,
34	! comment fields (!:) to be used for variable explanations added to
35	! all variable declaration statements
36	!
37	! 1318 2014-03-17 13:35:16Z raasch
38	! module interfaces removed
39	!
40	! 1306 2014-03-13 14:30:59Z raasch
41	! openmp sections removed from the overlap branch,
42	! second argument removed from parameter list
43	!
44	! 1216 2013-08-26 09:31:42Z raasch
45	! resorting of arrays moved to separate routines resort_for_...,
46	! one argument, used as temporary work array, removed from all transpose
47	! routines
48	! overlapping fft / transposition implemented
49	!
50	! 1212 2013-08-15 08:46:27Z raasch
51	! tridia routines moved to seperate module tridia_solver
52	!
53	! 1208 2013-08-13 06:41:49Z raasch
54	! acc-update clauses added for "ar" so that ffts other than cufft can also be
55	! used (although they are not ported and will give a poor performance)
56	!
57	! 1111 2013-03-08 23:54:10Z raasch
58	! further openACC porting of non-parallel (MPI) branch:
59	! tridiagonal routines split into extermal subroutines (instead using CONTAINS),
60	! no distinction between parallel/non-parallel in poisfft and tridia any more,
61	! tridia routines moved to end of file because of probable bug in PGI compiler 12.5
62	! (otherwise "invalid device function" is indicated during runtime),
63	! optimization of tridia routines: constant elements and coefficients of tri are
64	! stored in seperate arrays ddzuw and tric, last dimension of tri reduced from 5
65	! to 2,
66	! poisfft_init is now called internally from poisfft, maketri is called from
67	! poisfft_init,
68	! ibc_p_b = 2 removed
69	!
70	! 1106 2013-03-04 05:31:38Z raasch
71	! routines fftx, ffty, fftxp, fftyp removed, calls replaced by fft_x, fft_y,
72	! in the 1D-decomposition routines fft_x, ffty are replaced by fft_x_1d,
73	! fft_y_1d
74	!
75	! 1103 2013-02-20 02:15:53Z raasch
76	! tri, ar, and ar1 arguments in tridia-routines (2d) are removed because they
77	! sometimes cause segmentation faults with intel 12.1 compiler
78	!
79	! 1092 2013-02-02 11:24:22Z raasch
80	! unused variables removed
81	!
82	! 1036 2012-10-22 13:43:42Z raasch
83	! code put under GPL (PALM 3.9)
84	!
85	! 2012-09-21 07:03:55Z raasch
86	! FLOAT type conversion replaced by REAL
87	!
88	! 1003 2012-09-14 14:35:53Z raasch
89	! indices nxa, nya, etc. replaced by nx, ny, etc.
90	!
91	! 940 2012-07-09 14:31:00Z raasch
92	! special handling of tri-array as an argument in tridia_1dd routines switched
93	! off because it caused segmentation faults with intel 12.1 compiler
94	!
95	! 877 2012-04-03 11:21:44Z suehring
96	! Bugfix: Avoid divisions by zero in case of using a 'neumann' bc for the
97	! pressure at the top of the model domain.
98	!
99	! 809 2012-01-30 13:32:58Z maronga
100	! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
101	!
102	! 807 2012-01-25 11:53:51Z maronga
103	! New cpp directive "__check" implemented which is used by check_namelist_files
104	! (most of the code is unneeded by check_namelist_files).
105	!
106	! Revision 1.1 1997/07/24 11:24:14 raasch
107	! Initial revision
108	!
109	!
110	! Description:
111	! ------------
112	! Original version by Stephan Siano (pois3d), as of July 23, 1996
113	! Adapted for 2D-domain-decomposition by Siegfried Raasch, July 3, 1997
114	!
115	! Solves the Poisson equation with a 2D spectral method
116	! d^2 p / dx^2 + d^2 p / dy^2 + d^2 p / dz^2 = s
117	!
118	! Input:
119	! real ar contains (nnz,nny,nnx) elements of the velocity divergence,
120	! starting from (1,nys,nxl)
121	!
122	! Output:
123	! real ar contains the solution for perturbation pressure p
124	!------------------------------------------------------------------------------!
125
126	USE fft_xy, &
127	ONLY: fft_init, fft_y, fft_y_1d, fft_y_m, fft_x, fft_x_1d, fft_x_m
128
129	USE indices, &
130	ONLY: nnx, nny, nx, nxl, nxr, ny, nys, nyn, nz
131
132	USE transpose_indices, &
133	ONLY: nxl_y, nxl_z, nxr_y, nxr_z, nys_x, nys_z, nyn_x, nyn_z, nzb_x, &
134	nzb_y, nzt_x, nzt_y
135
136	USE tridia_solver, &
137	ONLY: tridia_1dd, tridia_init, tridia_substi, tridia_substi_overlap
138
139	IMPLICIT NONE
140
141	LOGICAL, SAVE :: poisfft_initialized = .FALSE.
142
143	PRIVATE
144
145	#if ! defined ( __check )
146	PUBLIC poisfft, poisfft_init
147
148	INTERFACE poisfft
149	MODULE PROCEDURE poisfft
150	END INTERFACE poisfft
151
152	INTERFACE poisfft_init
153	MODULE PROCEDURE poisfft_init
154	END INTERFACE poisfft_init
155	#else
156	PUBLIC poisfft_init
157
158	INTERFACE poisfft_init
159	MODULE PROCEDURE poisfft_init
160	END INTERFACE poisfft_init
161	#endif
162
163	CONTAINS
164
165	SUBROUTINE poisfft_init
166
167	USE arrays_3d, &
168	ONLY: ddzu_pres, ddzw
169
170	USE kinds
171
172	IMPLICIT NONE
173
174	INTEGER(iwp) :: k !:
175
176
177	CALL fft_init
178
179	CALL tridia_init
180
181	poisfft_initialized = .TRUE.
182
183	END SUBROUTINE poisfft_init
184
185
186	#if ! defined ( __check )
187	SUBROUTINE poisfft( ar )
188
189	USE control_parameters, &
190	ONLY: fft_method, transpose_compute_overlap
191
192	USE cpulog, &
193	ONLY: cpu_log, cpu_log_nowait, log_point_s
194
195	USE kinds
196
197	USE pegrid
198
199	IMPLICIT NONE
200
201	INTEGER(iwp) :: ii !:
202	INTEGER(iwp) :: iind !:
203	INTEGER(iwp) :: inew !:
204	INTEGER(iwp) :: jj !:
205	INTEGER(iwp) :: jind !:
206	INTEGER(iwp) :: jnew !:
207	INTEGER(iwp) :: ki !:
208	INTEGER(iwp) :: kk !:
209	INTEGER(iwp) :: knew !:
210	INTEGER(iwp) :: n !:
211	INTEGER(iwp) :: nblk !:
212	INTEGER(iwp) :: nnx_y !:
213	INTEGER(iwp) :: nny_z !:
214	INTEGER(iwp) :: nnz_t !:
215	INTEGER(iwp) :: nnz_x !:
216	INTEGER(iwp) :: nxl_y_bound !:
217	INTEGER(iwp) :: nxr_y_bound !:
218
219	INTEGER(iwp), DIMENSION(4) :: isave !:
220
221	REAL(wp), DIMENSION(1:nz,nys:nyn,nxl:nxr) :: ar !:
222	!$acc declare create( ar_inv )
223	REAL(wp), DIMENSION(nys:nyn,nxl:nxr,1:nz) :: ar_inv !:
224
225	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ar1 !:
226	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_in !:
227	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_inv !:
228	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_out_y !:
229	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f_out_z !:
230
231
232	CALL cpu_log( log_point_s(3), 'poisfft', 'start' )
233
234	IF ( .NOT. poisfft_initialized ) CALL poisfft_init
235
236	!
237	!-- Two-dimensional Fourier Transformation in x- and y-direction.
238	IF ( pdims(2) == 1 .AND. pdims(1) > 1 ) THEN
239
240	!
241	!-- 1d-domain-decomposition along x:
242	!-- FFT along y and transposition y --> x
243	CALL ffty_tr_yx( ar, ar )
244
245	!
246	!-- FFT along x, solving the tridiagonal system and backward FFT
247	CALL fftx_tri_fftx( ar )
248
249	!
250	!-- Transposition x --> y and backward FFT along y
251	CALL tr_xy_ffty( ar, ar )
252
253	ELSEIF ( pdims(1) == 1 .AND. pdims(2) > 1 ) THEN
254
255	!
256	!-- 1d-domain-decomposition along y:
257	!-- FFT along x and transposition x --> y
258	CALL fftx_tr_xy( ar, ar )
259
260	!
261	!-- FFT along y, solving the tridiagonal system and backward FFT
262	CALL ffty_tri_ffty( ar )
263
264	!
265	!-- Transposition y --> x and backward FFT along x
266	CALL tr_yx_fftx( ar, ar )
267
268	ELSEIF ( .NOT. transpose_compute_overlap ) THEN
269
270	!
271	!-- 2d-domain-decomposition or no decomposition (1 PE run)
272	!-- Transposition z --> x
273	CALL cpu_log( log_point_s(5), 'transpo forward', 'start' )
274	CALL resort_for_zx( ar, ar_inv )
275	CALL transpose_zx( ar_inv, ar )
276	CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
277
278	CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
279	IF ( fft_method /= 'system-specific' ) THEN
280	!$acc update host( ar )
281	ENDIF
282	CALL fft_x( ar, 'forward' )
283	IF ( fft_method /= 'system-specific' ) THEN
284	!$acc update device( ar )
285	ENDIF
286	CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
287
288	!
289	!-- Transposition x --> y
290	CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
291	CALL resort_for_xy( ar, ar_inv )
292	CALL transpose_xy( ar_inv, ar )
293	CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
294
295	CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
296	IF ( fft_method /= 'system-specific' ) THEN
297	!$acc update host( ar )
298	ENDIF
299	CALL fft_y( ar, 'forward', ar_tr = ar, &
300	nxl_y_bound = nxl_y, nxr_y_bound = nxr_y, &
301	nxl_y_l = nxl_y, nxr_y_l = nxr_y )
302	IF ( fft_method /= 'system-specific' ) THEN
303	!$acc update device( ar )
304	ENDIF
305	CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
306
307	!
308	!-- Transposition y --> z
309	CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
310	CALL resort_for_yz( ar, ar_inv )
311	CALL transpose_yz( ar_inv, ar )
312	CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
313
314	!
315	!-- Solve the tridiagonal equation system along z
316	CALL cpu_log( log_point_s(6), 'tridia', 'start' )
317	CALL tridia_substi( ar )
318	CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
319
320	!
321	!-- Inverse Fourier Transformation
322	!-- Transposition z --> y
323	CALL cpu_log( log_point_s(8), 'transpo invers', 'start' )
324	CALL transpose_zy( ar, ar_inv )
325	CALL resort_for_zy( ar_inv, ar )
326	CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
327
328	CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
329	IF ( fft_method /= 'system-specific' ) THEN
330	!$acc update host( ar )
331	ENDIF
332	CALL fft_y( ar, 'backward', ar_tr = ar, &
333	nxl_y_bound = nxl_y, nxr_y_bound = nxr_y, &
334	nxl_y_l = nxl_y, nxr_y_l = nxr_y )
335	IF ( fft_method /= 'system-specific' ) THEN
336	!$acc update device( ar )
337	ENDIF
338	CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
339
340	!
341	!-- Transposition y --> x
342	CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
343	CALL transpose_yx( ar, ar_inv )
344	CALL resort_for_yx( ar_inv, ar )
345	CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
346
347	CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
348	IF ( fft_method /= 'system-specific' ) THEN
349	!$acc update host( ar )
350	ENDIF
351	CALL fft_x( ar, 'backward' )
352	IF ( fft_method /= 'system-specific' ) THEN
353	!$acc update device( ar )
354	ENDIF
355	CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
356
357	!
358	!-- Transposition x --> z
359	CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
360	CALL transpose_xz( ar, ar_inv )
361	CALL resort_for_xz( ar_inv, ar )
362	CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
363
364	ELSE
365
366	!
367	!-- 2d-domain-decomposition or no decomposition (1 PE run) with
368	!-- overlapping transposition / fft
369	!-- cputime logging must not use barriers, which would prevent overlapping
370	ALLOCATE( f_out_y(0:ny,nxl_y:nxr_y,nzb_y:nzt_y), &
371	f_out_z(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
372	!
373	!-- Transposition z --> x + subsequent fft along x
374	ALLOCATE( f_inv(nys:nyn,nxl:nxr,1:nz) )
375	CALL resort_for_zx( ar, f_inv )
376	!
377	!-- Save original indices and gridpoint counter
378	isave(1) = nz
379	isave(2) = nzb_x
380	isave(3) = nzt_x
381	isave(4) = sendrecvcount_zx
382	!
383	!-- Set new indices for transformation
384	nblk = nz / pdims(1)
385	nz = pdims(1)
386	nnz_x = 1
387	nzb_x = 1 + myidx * nnz_x
388	nzt_x = ( myidx + 1 ) * nnz_x
389	sendrecvcount_zx = nnx * nny * nnz_x
390
391	ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
392	ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) )
393
394	DO kk = 1, nblk
395
396	IF ( kk == 1 ) THEN
397	CALL cpu_log( log_point_s(5), 'transpo forward', 'start', cpu_log_nowait )
398	ELSE
399	CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
400	ENDIF
401
402	DO knew = 1, nz
403	ki = kk + nblk * ( knew - 1 )
404	f_in(:,:,knew) = f_inv(:,:,ki)
405	ENDDO
406
407	CALL transpose_zx( f_in, ar1(:,:,:))
408	CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
409
410	IF ( kk == 1 ) THEN
411	CALL cpu_log( log_point_s(4), 'fft_x', 'start', cpu_log_nowait )
412	ELSE
413	CALL cpu_log( log_point_s(4), 'fft_x', 'continue', cpu_log_nowait )
414	ENDIF
415
416	n = isave(2) + kk - 1
417	CALL fft_x( ar1(:,:,:), 'forward', ar_2d = f_out_z(:,:,n))
418	CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
419
420	ENDDO
421	!
422	!-- Restore original indices/counters
423	nz = isave(1)
424	nzb_x = isave(2)
425	nzt_x = isave(3)
426	sendrecvcount_zx = isave(4)
427
428	DEALLOCATE( ar1, f_in, f_inv )
429
430	!
431	!-- Transposition x --> y + subsequent fft along y
432	ALLOCATE( f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
433	CALL resort_for_xy( f_out_z, f_inv )
434	!
435	!-- Save original indices and gridpoint counter
436	isave(1) = nx
437	isave(2) = nxl_y
438	isave(3) = nxr_y
439	isave(4) = sendrecvcount_xy
440	!
441	!-- Set new indices for transformation
442	nblk = ( ( nx+1 ) / pdims(2) ) - 1
443	nx = pdims(2)
444	nnx_y = 1
445	nxl_y = myidy * nnx_y
446	nxr_y = ( myidy + 1 ) * nnx_y - 1
447	sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 )
448
449	ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) )
450	ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
451
452	DO ii = 0, nblk
453
454	CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
455
456	DO inew = 0, nx-1
457	iind = ii + ( nblk + 1 ) * inew
458	f_in(:,:,inew) = f_inv(:,:,iind)
459	ENDDO
460
461	CALL transpose_xy( f_in, ar1(:,:,:) )
462
463	CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
464
465	IF ( ii == 1 ) THEN
466	CALL cpu_log( log_point_s(7), 'fft_y', 'start', cpu_log_nowait )
467	ELSE
468	CALL cpu_log( log_point_s(7), 'fft_y', 'continue', cpu_log_nowait )
469	ENDIF
470
471	nxl_y_bound = isave(2)
472	nxr_y_bound = isave(3)
473	n = isave(2) + ii
474	CALL fft_y( ar1(:,:,:), 'forward', ar_tr = f_out_y, &
475	nxl_y_bound = nxl_y_bound, nxr_y_bound = nxr_y_bound, &
476	nxl_y_l = n, nxr_y_l = n )
477
478	CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
479
480	ENDDO
481	!
482	!-- Restore original indices/counters
483	nx = isave(1)
484	nxl_y = isave(2)
485	nxr_y = isave(3)
486	sendrecvcount_xy = isave(4)
487
488	DEALLOCATE( ar1, f_in, f_inv )
489
490	!
491	!-- Transposition y --> z + subsequent tridia + resort for z --> y
492	ALLOCATE( f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) )
493	CALL resort_for_yz( f_out_y, f_inv )
494	!
495	!-- Save original indices and gridpoint counter
496	isave(1) = ny
497	isave(2) = nys_z
498	isave(3) = nyn_z
499	isave(4) = sendrecvcount_yz
500	!
501	!-- Set new indices for transformation
502	nblk = ( ( ny+1 ) / pdims(1) ) - 1
503	ny = pdims(1)
504	nny_z = 1
505	nys_z = myidx * nny_z
506	nyn_z = ( myidx + 1 ) * nny_z - 1
507	sendrecvcount_yz = ( nxr_y-nxl_y+1 ) * nny_z * ( nzt_y-nzb_y+1 )
508
509	ALLOCATE( ar1(nxl_z:nxr_z,nys_z:nyn_z,1:nz) )
510	ALLOCATE( f_in(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) )
511
512	DO jj = 0, nblk
513	!
514	!-- Forward Fourier Transformation
515	!-- Transposition y --> z
516	CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
517
518	DO jnew = 0, ny-1
519	jind = jj + ( nblk + 1 ) * jnew
520	f_in(:,:,jnew) = f_inv(:,:,jind)
521	ENDDO
522
523	CALL transpose_yz( f_in, ar1(:,:,:) )
524
525	IF ( jj == nblk ) THEN
526	CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
527	ELSE
528	CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
529	ENDIF
530
531	!
532	!-- Solve the tridiagonal equation system along z
533	CALL cpu_log( log_point_s(6), 'tridia', 'start', cpu_log_nowait )
534
535	n = isave(2) + jj
536	CALL tridia_substi_overlap( ar1(:,:,:), n )
537
538	CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
539
540	!
541	!-- Inverse Fourier Transformation
542	!-- Transposition z --> y
543	!-- Only one thread should call MPI routines, therefore forward and
544	!-- backward tranpose are in the same section
545	IF ( jj == 0 ) THEN
546	CALL cpu_log( log_point_s(8), 'transpo invers', 'start', cpu_log_nowait )
547	ELSE
548	CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
549	ENDIF
550
551	CALL transpose_zy( ar1(:,:,:), f_in )
552
553	DO jnew = 0, ny-1
554	jind = jj + ( nblk + 1 ) * jnew
555	f_inv(:,:,jind) = f_in(:,:,jnew)
556	ENDDO
557
558	CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
559
560	ENDDO
561	!
562	!-- Restore original indices/counters
563	ny = isave(1)
564	nys_z = isave(2)
565	nyn_z = isave(3)
566	sendrecvcount_yz = isave(4)
567
568	CALL resort_for_zy( f_inv, f_out_y )
569
570	DEALLOCATE( ar1, f_in, f_inv )
571
572	!
573	!-- fft along y backward + subsequent transposition y --> x
574	ALLOCATE( f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
575	!
576	!-- Save original indices and gridpoint counter
577	isave(1) = nx
578	isave(2) = nxl_y
579	isave(3) = nxr_y
580	isave(4) = sendrecvcount_xy
581	!
582	!-- Set new indices for transformation
583	nblk = (( nx+1 ) / pdims(2) ) - 1
584	nx = pdims(2)
585	nnx_y = 1
586	nxl_y = myidy * nnx_y
587	nxr_y = ( myidy + 1 ) * nnx_y - 1
588	sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 )
589
590	ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) )
591	ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
592
593	DO ii = 0, nblk
594
595	CALL cpu_log( log_point_s(7), 'fft_y', 'continue', cpu_log_nowait )
596
597	n = isave(2) + ii
598	nxl_y_bound = isave(2)
599	nxr_y_bound = isave(3)
600
601	CALL fft_y( ar1(:,:,:), 'backward', ar_tr = f_out_y, &
602	nxl_y_bound = nxl_y_bound, nxr_y_bound = nxr_y_bound, &
603	nxl_y_l = n, nxr_y_l = n )
604
605	IF ( ii == nblk ) THEN
606	CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
607	ELSE
608	CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
609	ENDIF
610
611	CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
612
613	CALL transpose_yx( ar1(:,:,:), f_in )
614
615	DO inew = 0, nx-1
616	iind = ii + (nblk+1) * inew
617	f_inv(:,:,iind) = f_in(:,:,inew)
618	ENDDO
619
620	CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
621
622	ENDDO
623	!
624	!-- Restore original indices/counters
625	nx = isave(1)
626	nxl_y = isave(2)
627	nxr_y = isave(3)
628	sendrecvcount_xy = isave(4)
629
630	CALL resort_for_yx( f_inv, f_out_z )
631
632	DEALLOCATE( ar1, f_in, f_inv )
633
634	!
635	!-- fft along x backward + subsequent final transposition x --> z
636	ALLOCATE( f_inv(nys:nyn,nxl:nxr,1:nz) )
637	!
638	!-- Save original indices and gridpoint counter
639	isave(1) = nz
640	isave(2) = nzb_x
641	isave(3) = nzt_x
642	isave(4) = sendrecvcount_zx
643	!
644	!-- Set new indices for transformation
645	nblk = nz / pdims(1)
646	nz = pdims(1)
647	nnz_x = 1
648	nzb_x = 1 + myidx * nnz_x
649	nzt_x = ( myidx + 1 ) * nnz_x
650	sendrecvcount_zx = nnx * nny * nnz_x
651
652	ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
653	ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) )
654
655	DO kk = 1, nblk
656
657	CALL cpu_log( log_point_s(4), 'fft_x', 'continue', cpu_log_nowait )
658
659	n = isave(2) + kk - 1
660	CALL fft_x( ar1(:,:,:), 'backward', f_out_z(:,:,n))
661
662	IF ( kk == nblk ) THEN
663	CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
664	ELSE
665	CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
666	ENDIF
667
668	CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
669
670	CALL transpose_xz( ar1(:,:,:), f_in )
671
672	DO knew = 1, nz
673	ki = kk + nblk * (knew-1)
674	f_inv(:,:,ki) = f_in(:,:,knew)
675	ENDDO
676
677	IF ( kk == nblk ) THEN
678	CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
679	ELSE
680	CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
681	ENDIF
682
683	ENDDO
684	!
685	!-- Restore original indices/counters
686	nz = isave(1)
687	nzb_x = isave(2)
688	nzt_x = isave(3)
689	sendrecvcount_zx = isave(4)
690
691	CALL resort_for_xz( f_inv, ar )
692
693	DEALLOCATE( ar1, f_in, f_inv )
694
695	ENDIF
696
697	CALL cpu_log( log_point_s(3), 'poisfft', 'stop' )
698
699	END SUBROUTINE poisfft
700
701
702
703	SUBROUTINE ffty_tr_yx( f_in, f_out )
704
705	!------------------------------------------------------------------------------!
706	! Fourier-transformation along y with subsequent transposition y --> x for
707	! a 1d-decomposition along x
708	!
709	! ATTENTION: The performance of this routine is much faster on the NEC-SX6,
710	! if the first index of work_ffty_vec is odd. Otherwise
711	! memory bank conflicts may occur (especially if the index is a
712	! multiple of 128). That's why work_ffty_vec is dimensioned as
713	! 0:ny+1.
714	! Of course, this will not work if users are using an odd number
715	! of gridpoints along y.
716	!------------------------------------------------------------------------------!
717
718	USE control_parameters, &
719	ONLY: host
720
721	USE cpulog, &
722	ONLY: cpu_log, log_point_s
723
724	USE kinds
725
726	USE pegrid
727
728	IMPLICIT NONE
729
730	INTEGER(iwp) :: i !:
731	INTEGER(iwp) :: iend !:
732	INTEGER(iwp) :: iouter !:
733	INTEGER(iwp) :: ir !:
734	INTEGER(iwp) :: j !:
735	INTEGER(iwp) :: k !:
736
737	INTEGER(iwp), PARAMETER :: stridex = 4 !:
738
739	REAL(wp), DIMENSION(0:ny,stridex) :: work_ffty !:
740	#if defined( __nec )
741	REAL(wp), DIMENSION(0:ny+1,1:nz,nxl:nxr) :: work_ffty_vec !:
742	#endif
743	REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr) :: f_in !:
744	REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) :: f_out !:
745	REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny) :: work !:
746
747	!
748	!-- Carry out the FFT along y, where all data are present due to the
749	!-- 1d-decomposition along x. Resort the data in a way that x becomes
750	!-- the first index.
751	CALL cpu_log( log_point_s(7), 'fft_y_1d', 'start' )
752
753	IF ( host(1:3) == 'nec' ) THEN
754	#if defined( __nec )
755	!
756	!-- Code optimized for vector processors
757	!$OMP PARALLEL PRIVATE ( i, j, k )
758	!$OMP DO
759	DO i = nxl, nxr
760
761	DO j = 0, ny
762	DO k = 1, nz
763	work_ffty_vec(j,k,i) = f_in(k,j,i)
764	ENDDO
765	ENDDO
766
767	CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'forward' )
768
769	ENDDO
770
771	!$OMP DO
772	DO k = 1, nz
773	DO j = 0, ny
774	DO i = nxl, nxr
775	work(i,k,j) = work_ffty_vec(j,k,i)
776	ENDDO
777	ENDDO
778	ENDDO
779	!$OMP END PARALLEL
780	#endif
781
782	ELSE
783
784	!
785	!-- Cache optimized code.
786	!-- The i-(x-)direction is split into a strided outer loop and an inner
787	!-- loop for better cache performance
788	!$OMP PARALLEL PRIVATE (i,iend,iouter,ir,j,k,work_ffty)
789	!$OMP DO
790	DO iouter = nxl, nxr, stridex
791
792	iend = MIN( iouter+stridex-1, nxr ) ! Upper bound for inner i loop
793
794	DO k = 1, nz
795
796	DO i = iouter, iend
797
798	ir = i-iouter+1 ! counter within a stride
799	DO j = 0, ny
800	work_ffty(j,ir) = f_in(k,j,i)
801	ENDDO
802	!
803	!-- FFT along y
804	CALL fft_y_1d( work_ffty(:,ir), 'forward' )
805
806	ENDDO
807
808	!
809	!-- Resort
810	DO j = 0, ny
811	DO i = iouter, iend
812	work(i,k,j) = work_ffty(j,i-iouter+1)
813	ENDDO
814	ENDDO
815
816	ENDDO
817
818	ENDDO
819	!$OMP END PARALLEL
820
821	ENDIF
822	CALL cpu_log( log_point_s(7), 'fft_y_1d', 'pause' )
823
824	!
825	!-- Transpose array
826	#if defined( __parallel )
827	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
828	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
829	CALL MPI_ALLTOALL( work(nxl,1,0), sendrecvcount_xy, MPI_REAL, &
830	f_out(1,1,nys_x,1), sendrecvcount_xy, MPI_REAL, &
831	comm1dx, ierr )
832	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
833	#endif
834
835	END SUBROUTINE ffty_tr_yx
836
837
838	SUBROUTINE tr_xy_ffty( f_in, f_out )
839
840	!------------------------------------------------------------------------------!
841	! Transposition x --> y with a subsequent backward Fourier transformation for
842	! a 1d-decomposition along x
843	!------------------------------------------------------------------------------!
844
845	USE control_parameters, &
846	ONLY: host
847
848	USE cpulog, &
849	ONLY: cpu_log, log_point_s
850
851	USE kinds
852
853	USE pegrid
854
855	IMPLICIT NONE
856
857	INTEGER(iwp) :: i !:
858	INTEGER(iwp) :: iend !:
859	INTEGER(iwp) :: iouter !:
860	INTEGER(iwp) :: ir !:
861	INTEGER(iwp) :: j !:
862	INTEGER(iwp) :: k !:
863
864	INTEGER(iwp), PARAMETER :: stridex = 4 !:
865
866	REAL(wp), DIMENSION(0:ny,stridex) :: work_ffty !:
867	#if defined( __nec )
868	REAL(wp), DIMENSION(0:ny+1,1:nz,nxl:nxr) :: work_ffty_vec !:
869	#endif
870	REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) :: f_in !:
871	REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr) :: f_out !:
872	REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny) :: work !:
873
874	!
875	!-- Transpose array
876	#if defined( __parallel )
877	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
878	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
879	CALL MPI_ALLTOALL( f_in(1,1,nys_x,1), sendrecvcount_xy, MPI_REAL, &
880	work(nxl,1,0), sendrecvcount_xy, MPI_REAL, &
881	comm1dx, ierr )
882	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
883	#endif
884
885	!
886	!-- Resort the data in a way that y becomes the first index and carry out the
887	!-- backward fft along y.
888	CALL cpu_log( log_point_s(7), 'fft_y_1d', 'continue' )
889
890	IF ( host(1:3) == 'nec' ) THEN
891	#if defined( __nec )
892	!
893	!-- Code optimized for vector processors
894	!$OMP PARALLEL PRIVATE ( i, j, k )
895	!$OMP DO
896	DO k = 1, nz
897	DO j = 0, ny
898	DO i = nxl, nxr
899	work_ffty_vec(j,k,i) = work(i,k,j)
900	ENDDO
901	ENDDO
902	ENDDO
903
904	!$OMP DO
905	DO i = nxl, nxr
906
907	CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'backward' )
908
909	DO j = 0, ny
910	DO k = 1, nz
911	f_out(k,j,i) = work_ffty_vec(j,k,i)
912	ENDDO
913	ENDDO
914
915	ENDDO
916	!$OMP END PARALLEL
917	#endif
918
919	ELSE
920
921	!
922	!-- Cache optimized code.
923	!-- The i-(x-)direction is split into a strided outer loop and an inner
924	!-- loop for better cache performance
925	!$OMP PARALLEL PRIVATE ( i, iend, iouter, ir, j, k, work_ffty )
926	!$OMP DO
927	DO iouter = nxl, nxr, stridex
928
929	iend = MIN( iouter+stridex-1, nxr ) ! Upper bound for inner i loop
930
931	DO k = 1, nz
932	!
933	!-- Resort
934	DO j = 0, ny
935	DO i = iouter, iend
936	work_ffty(j,i-iouter+1) = work(i,k,j)
937	ENDDO
938	ENDDO
939
940	DO i = iouter, iend
941
942	!
943	!-- FFT along y
944	ir = i-iouter+1 ! counter within a stride
945	CALL fft_y_1d( work_ffty(:,ir), 'backward' )
946
947	DO j = 0, ny
948	f_out(k,j,i) = work_ffty(j,ir)
949	ENDDO
950	ENDDO
951
952	ENDDO
953
954	ENDDO
955	!$OMP END PARALLEL
956
957	ENDIF
958
959	CALL cpu_log( log_point_s(7), 'fft_y_1d', 'stop' )
960
961	END SUBROUTINE tr_xy_ffty
962
963
964	SUBROUTINE fftx_tri_fftx( ar )
965
966	!------------------------------------------------------------------------------!
967	! FFT along x, solution of the tridiagonal system and backward FFT for
968	! a 1d-decomposition along x
969	!
970	! WARNING: this subroutine may still not work for hybrid parallelization
971	! with OpenMP (for possible necessary changes see the original
972	! routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002)
973	!------------------------------------------------------------------------------!
974
975	USE control_parameters, &
976	ONLY: host
977
978	USE cpulog, &
979	ONLY: cpu_log, log_point_s
980
981	USE grid_variables, &
982	ONLY: ddx2, ddy2
983
984	USE kinds
985
986	USE pegrid
987
988	IMPLICIT NONE
989
990	INTEGER(iwp) :: i !:
991	INTEGER(iwp) :: j !:
992	INTEGER(iwp) :: k !:
993	INTEGER(iwp) :: m !:
994	INTEGER(iwp) :: n !:
995	INTEGER(iwp) :: omp_get_thread_num !:
996	INTEGER(iwp) :: tn !:
997
998	REAL(wp), DIMENSION(0:nx) :: work_fftx !:
999	REAL(wp), DIMENSION(0:nx,1:nz) :: work_trix !:
1000	REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) :: ar !:
1001	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tri !:
1002
1003
1004	CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'start' )
1005
1006	ALLOCATE( tri(5,0:nx,0:nz-1,0:threads_per_task-1) )
1007
1008	tn = 0 ! Default thread number in case of one thread
1009	!$OMP PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_fftx, work_trix )
1010	DO j = nys_x, nyn_x
1011
1012	!$ tn = omp_get_thread_num()
1013
1014	IF ( host(1:3) == 'nec' ) THEN
1015	!
1016	!-- Code optimized for vector processors
1017	DO k = 1, nz
1018
1019	m = 0
1020	DO n = 1, pdims(1)
1021	DO i = 1, nnx
1022	work_trix(m,k) = ar(i,k,j,n)
1023	m = m + 1
1024	ENDDO
1025	ENDDO
1026
1027	ENDDO
1028
1029	CALL fft_x_m( work_trix, 'forward' )
1030
1031	ELSE
1032	!
1033	!-- Cache optimized code
1034	DO k = 1, nz
1035
1036	m = 0
1037	DO n = 1, pdims(1)
1038	DO i = 1, nnx
1039	work_fftx(m) = ar(i,k,j,n)
1040	m = m + 1
1041	ENDDO
1042	ENDDO
1043
1044	CALL fft_x_1d( work_fftx, 'forward' )
1045
1046	DO i = 0, nx
1047	work_trix(i,k) = work_fftx(i)
1048	ENDDO
1049
1050	ENDDO
1051
1052	ENDIF
1053
1054	!
1055	!-- Solve the linear equation system
1056	CALL tridia_1dd( ddx2, ddy2, nx, ny, j, work_trix, tri(:,:,:,tn) )
1057
1058	IF ( host(1:3) == 'nec' ) THEN
1059	!
1060	!-- Code optimized for vector processors
1061	CALL fft_x_m( work_trix, 'backward' )
1062
1063	DO k = 1, nz
1064
1065	m = 0
1066	DO n = 1, pdims(1)
1067	DO i = 1, nnx
1068	ar(i,k,j,n) = work_trix(m,k)
1069	m = m + 1
1070	ENDDO
1071	ENDDO
1072
1073	ENDDO
1074
1075	ELSE
1076	!
1077	!-- Cache optimized code
1078	DO k = 1, nz
1079
1080	DO i = 0, nx
1081	work_fftx(i) = work_trix(i,k)
1082	ENDDO
1083
1084	CALL fft_x_1d( work_fftx, 'backward' )
1085
1086	m = 0
1087	DO n = 1, pdims(1)
1088	DO i = 1, nnx
1089	ar(i,k,j,n) = work_fftx(m)
1090	m = m + 1
1091	ENDDO
1092	ENDDO
1093
1094	ENDDO
1095
1096	ENDIF
1097
1098	ENDDO
1099
1100	DEALLOCATE( tri )
1101
1102	CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'stop' )
1103
1104	END SUBROUTINE fftx_tri_fftx
1105
1106
1107	SUBROUTINE fftx_tr_xy( f_in, f_out )
1108
1109	!------------------------------------------------------------------------------!
1110	! Fourier-transformation along x with subsequent transposition x --> y for
1111	! a 1d-decomposition along y
1112	!
1113	! ATTENTION: The NEC-branch of this routine may significantly profit from
1114	! further optimizations. So far, performance is much worse than
1115	! for routine ffty_tr_yx (more than three times slower).
1116	!------------------------------------------------------------------------------!
1117
1118	USE control_parameters, &
1119	ONLY: host
1120
1121	USE cpulog, &
1122	ONLY: cpu_log, log_point_s
1123
1124	USE kinds
1125
1126	USE pegrid
1127
1128	IMPLICIT NONE
1129
1130	INTEGER(iwp) :: i !:
1131	INTEGER(iwp) :: j !:
1132	INTEGER(iwp) :: k !:
1133
1134	REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn) :: work_fftx !:
1135	REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx) :: f_in !:
1136	REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) :: f_out !:
1137	REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx) :: work !:
1138
1139	!
1140	!-- Carry out the FFT along x, where all data are present due to the
1141	!-- 1d-decomposition along y. Resort the data in a way that y becomes
1142	!-- the first index.
1143	CALL cpu_log( log_point_s(4), 'fft_x_1d', 'start' )
1144
1145	IF ( host(1:3) == 'nec' ) THEN
1146	!
1147	!-- Code for vector processors
1148	!$OMP PARALLEL PRIVATE ( i, j, k )
1149	!$OMP DO
1150	DO i = 0, nx
1151
1152	DO j = nys, nyn
1153	DO k = 1, nz
1154	work_fftx(i,k,j) = f_in(k,j,i)
1155	ENDDO
1156	ENDDO
1157
1158	ENDDO
1159
1160	!$OMP DO
1161	DO j = nys, nyn
1162
1163	CALL fft_x_m( work_fftx(:,:,j), 'forward' )
1164
1165	DO k = 1, nz
1166	DO i = 0, nx
1167	work(j,k,i) = work_fftx(i,k,j)
1168	ENDDO
1169	ENDDO
1170
1171	ENDDO
1172	!$OMP END PARALLEL
1173
1174	ELSE
1175
1176	!
1177	!-- Cache optimized code (there might be still a potential for better
1178	!-- optimization).
1179	!$OMP PARALLEL PRIVATE (i,j,k)
1180	!$OMP DO
1181	DO i = 0, nx
1182
1183	DO j = nys, nyn
1184	DO k = 1, nz
1185	work_fftx(i,k,j) = f_in(k,j,i)
1186	ENDDO
1187	ENDDO
1188
1189	ENDDO
1190
1191	!$OMP DO
1192	DO j = nys, nyn
1193	DO k = 1, nz
1194
1195	CALL fft_x_1d( work_fftx(0:nx,k,j), 'forward' )
1196
1197	DO i = 0, nx
1198	work(j,k,i) = work_fftx(i,k,j)
1199	ENDDO
1200	ENDDO
1201
1202	ENDDO
1203	!$OMP END PARALLEL
1204
1205	ENDIF
1206	CALL cpu_log( log_point_s(4), 'fft_x_1d', 'pause' )
1207
1208	!
1209	!-- Transpose array
1210	#if defined( __parallel )
1211	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1212	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1213	CALL MPI_ALLTOALL( work(nys,1,0), sendrecvcount_xy, MPI_REAL, &
1214	f_out(1,1,nxl_y,1), sendrecvcount_xy, MPI_REAL, &
1215	comm1dy, ierr )
1216	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1217	#endif
1218
1219	END SUBROUTINE fftx_tr_xy
1220
1221
1222	SUBROUTINE tr_yx_fftx( f_in, f_out )
1223
1224	!------------------------------------------------------------------------------!
1225	! Transposition y --> x with a subsequent backward Fourier transformation for
1226	! a 1d-decomposition along x
1227	!------------------------------------------------------------------------------!
1228
1229	USE control_parameters, &
1230	ONLY: host
1231
1232	USE cpulog, &
1233	ONLY: cpu_log, log_point_s
1234
1235	USE kinds
1236
1237	USE pegrid
1238
1239	IMPLICIT NONE
1240
1241	INTEGER(iwp) :: i !:
1242	INTEGER(iwp) :: j !:
1243	INTEGER(iwp) :: k !:
1244
1245	REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn) :: work_fftx !:
1246	REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) :: f_in !:
1247	REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx) :: f_out !:
1248	REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx) :: work !:
1249
1250	!
1251	!-- Transpose array
1252	#if defined( __parallel )
1253	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1254	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1255	CALL MPI_ALLTOALL( f_in(1,1,nxl_y,1), sendrecvcount_xy, MPI_REAL, &
1256	work(nys,1,0), sendrecvcount_xy, MPI_REAL, &
1257	comm1dy, ierr )
1258	CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1259	#endif
1260
1261	!
1262	!-- Carry out the FFT along x, where all data are present due to the
1263	!-- 1d-decomposition along y. Resort the data in a way that y becomes
1264	!-- the first index.
1265	CALL cpu_log( log_point_s(4), 'fft_x_1d', 'continue' )
1266
1267	IF ( host(1:3) == 'nec' ) THEN
1268	!
1269	!-- Code optimized for vector processors
1270	!$OMP PARALLEL PRIVATE ( i, j, k )
1271	!$OMP DO
1272	DO j = nys, nyn
1273
1274	DO k = 1, nz
1275	DO i = 0, nx
1276	work_fftx(i,k,j) = work(j,k,i)
1277	ENDDO
1278	ENDDO
1279
1280	CALL fft_x_m( work_fftx(:,:,j), 'backward' )
1281
1282	ENDDO
1283
1284	!$OMP DO
1285	DO i = 0, nx
1286	DO j = nys, nyn
1287	DO k = 1, nz
1288	f_out(k,j,i) = work_fftx(i,k,j)
1289	ENDDO
1290	ENDDO
1291	ENDDO
1292	!$OMP END PARALLEL
1293
1294	ELSE
1295
1296	!
1297	!-- Cache optimized code (there might be still a potential for better
1298	!-- optimization).
1299	!$OMP PARALLEL PRIVATE (i,j,k)
1300	!$OMP DO
1301	DO j = nys, nyn
1302	DO k = 1, nz
1303
1304	DO i = 0, nx
1305	work_fftx(i,k,j) = work(j,k,i)
1306	ENDDO
1307
1308	CALL fft_x_1d( work_fftx(0:nx,k,j), 'backward' )
1309
1310	ENDDO
1311	ENDDO
1312
1313	!$OMP DO
1314	DO i = 0, nx
1315	DO j = nys, nyn
1316	DO k = 1, nz
1317	f_out(k,j,i) = work_fftx(i,k,j)
1318	ENDDO
1319	ENDDO
1320	ENDDO
1321	!$OMP END PARALLEL
1322
1323	ENDIF
1324	CALL cpu_log( log_point_s(4), 'fft_x_1d', 'stop' )
1325
1326	END SUBROUTINE tr_yx_fftx
1327
1328
1329	SUBROUTINE ffty_tri_ffty( ar )
1330
1331	!------------------------------------------------------------------------------!
1332	! FFT along y, solution of the tridiagonal system and backward FFT for
1333	! a 1d-decomposition along y
1334	!
1335	! WARNING: this subroutine may still not work for hybrid parallelization
1336	! with OpenMP (for possible necessary changes see the original
1337	! routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002)
1338	!------------------------------------------------------------------------------!
1339
1340	USE control_parameters, &
1341	ONLY: host
1342
1343	USE cpulog, &
1344	ONLY: cpu_log, log_point_s
1345
1346	USE grid_variables, &
1347	ONLY: ddx2, ddy2
1348
1349	USE kinds
1350
1351	USE pegrid
1352
1353	IMPLICIT NONE
1354
1355	INTEGER(iwp) :: i !:
1356	INTEGER(iwp) :: j !:
1357	INTEGER(iwp) :: k !:
1358	INTEGER(iwp) :: m !:
1359	INTEGER(iwp) :: n !:
1360	INTEGER(iwp) :: omp_get_thread_num !:
1361	INTEGER(iwp) :: tn !:
1362
1363	REAL(wp), DIMENSION(0:ny) :: work_ffty !:
1364	REAL(wp), DIMENSION(0:ny,1:nz) :: work_triy !:
1365	REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) :: ar !:
1366	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tri !:
1367
1368
1369	CALL cpu_log( log_point_s(39), 'fft_y_1d + tridia', 'start' )
1370
1371	ALLOCATE( tri(5,0:ny,0:nz-1,0:threads_per_task-1) )
1372
1373	tn = 0 ! Default thread number in case of one thread
1374	!$OMP PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_ffty, work_triy )
1375	DO i = nxl_y, nxr_y
1376
1377	!$ tn = omp_get_thread_num()
1378
1379	IF ( host(1:3) == 'nec' ) THEN
1380	!
1381	!-- Code optimized for vector processors
1382	DO k = 1, nz
1383
1384	m = 0
1385	DO n = 1, pdims(2)
1386	DO j = 1, nny
1387	work_triy(m,k) = ar(j,k,i,n)
1388	m = m + 1
1389	ENDDO
1390	ENDDO
1391
1392	ENDDO
1393
1394	CALL fft_y_m( work_triy, ny, 'forward' )
1395
1396	ELSE
1397	!
1398	!-- Cache optimized code
1399	DO k = 1, nz
1400
1401	m = 0
1402	DO n = 1, pdims(2)
1403	DO j = 1, nny
1404	work_ffty(m) = ar(j,k,i,n)
1405	m = m + 1
1406	ENDDO
1407	ENDDO
1408
1409	CALL fft_y_1d( work_ffty, 'forward' )
1410
1411	DO j = 0, ny
1412	work_triy(j,k) = work_ffty(j)
1413	ENDDO
1414
1415	ENDDO
1416
1417	ENDIF
1418
1419	!
1420	!-- Solve the linear equation system
1421	CALL tridia_1dd( ddy2, ddx2, ny, nx, i, work_triy, tri(:,:,:,tn) )
1422
1423	IF ( host(1:3) == 'nec' ) THEN
1424	!
1425	!-- Code optimized for vector processors
1426	CALL fft_y_m( work_triy, ny, 'backward' )
1427
1428	DO k = 1, nz
1429
1430	m = 0
1431	DO n = 1, pdims(2)
1432	DO j = 1, nny
1433	ar(j,k,i,n) = work_triy(m,k)
1434	m = m + 1
1435	ENDDO
1436	ENDDO
1437
1438	ENDDO
1439
1440	ELSE
1441	!
1442	!-- Cache optimized code
1443	DO k = 1, nz
1444
1445	DO j = 0, ny
1446	work_ffty(j) = work_triy(j,k)
1447	ENDDO
1448
1449	CALL fft_y_1d( work_ffty, 'backward' )
1450
1451	m = 0
1452	DO n = 1, pdims(2)
1453	DO j = 1, nny
1454	ar(j,k,i,n) = work_ffty(m)
1455	m = m + 1
1456	ENDDO
1457	ENDDO
1458
1459	ENDDO
1460
1461	ENDIF
1462
1463	ENDDO
1464
1465	DEALLOCATE( tri )
1466
1467	CALL cpu_log( log_point_s(39), 'fft_y_1d + tridia', 'stop' )
1468
1469	END SUBROUTINE ffty_tri_ffty
1470
1471	#endif
1472
1473	END MODULE poisfft_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |