Home

Context Navigation

source: palm/trunk/SOURCE/fft_xy_mod.f90 @ 4646

Last change on this file since 4646 was 4646, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 59.0 KB

Line
1	!> @file fft_xy_mod.f90
2	!--------------------------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
7	! (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11	! Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
14	! <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2020 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------------------------!
18	!
19	! Current revisions:
20	! -----------------
21	!
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: fft_xy_mod.f90 4646 2020-08-24 16:02:40Z raasch $
26	! file re-formatted to follow the PALM coding standard
27	!
28	! 4370 2020-01-10 14:00:44Z raasch
29	! bugfix for Temperton-fft usage on GPU
30	!
31	! 4366 2020-01-09 08:12:43Z raasch
32	! Vectorized Temperton-fft added
33	!
34	! 4360 2020-01-07 11:25:50Z suehring
35	! Corrected "Former revisions" section
36	!
37	! 4069 2019-07-01 14:05:51Z Giersch
38	! Code added to avoid compiler warnings
39	!
40	! 3655 2019-01-07 16:51:22Z knoop
41	! OpenACC port for SPEC
42	!
43	! Revision 1.1 2002/06/11 13:00:49 raasch
44	! Initial revision
45	!
46	!
47	! Description:
48	! ------------
49	!> Fast Fourier transformation along x and y for 1d domain decomposition along x.
50	!> Original version: Klaus Ketelsen (May 2002)
51	!> @todo openmp support for vectorized Temperton fft
52	!------------------------------------------------------------------------------!
53	MODULE fft_xy
54
55
56	USE control_parameters, &
57	ONLY: fft_method, loop_optimization, message_string
58
59	USE cuda_fft_interfaces
60
61	USE indices, &
62	ONLY: nx, ny, nz
63
64	#if defined( __cuda_fft )
65	USE ISO_C_BINDING
66	#elif defined( __fftw )
67	USE, INTRINSIC :: ISO_C_BINDING
68	#endif
69
70	USE kinds
71
72	USE singleton, &
73	ONLY: fftn
74
75	USE temperton_fft
76
77	USE transpose_indices, &
78	ONLY: nxl_y, nxr_y, nyn_x, nys_x, nzb_x, nzb_y, nzt_x, nzt_y
79
80	IMPLICIT NONE
81
82	PRIVATE
83	PUBLIC fft_init, f_vec_x, fft_x, fft_x_1d, fft_x_m, fft_y, fft_y_1d, fft_y_m, temperton_fft_vec
84
85	INTEGER(iwp), DIMENSION(:), ALLOCATABLE, SAVE :: ifax_x !<
86	INTEGER(iwp), DIMENSION(:), ALLOCATABLE, SAVE :: ifax_y !<
87
88	LOGICAL, SAVE :: init_fft = .FALSE. !<
89	LOGICAL, SAVE :: temperton_fft_vec = .FALSE. !<
90
91	REAL(wp), SAVE :: dnx !<
92	REAL(wp), SAVE :: dny !<
93	REAL(wp), SAVE :: sqr_dnx !<
94	REAL(wp), SAVE :: sqr_dny !<
95
96	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: trigs_x !<
97	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: trigs_y !<
98
99	REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: f_vec_x
100
101	#if defined( __ibm )
102	INTEGER(iwp), PARAMETER :: nau1 = 20000 !<
103	INTEGER(iwp), PARAMETER :: nau2 = 22000 !<
104	!
105	!-- The following working arrays contain tables and have to be "save" and shared in OpenMP sense
106	REAL(wp), DIMENSION(nau1), SAVE :: aux1 !<
107	REAL(wp), DIMENSION(nau1), SAVE :: auy1 !<
108	REAL(wp), DIMENSION(nau1), SAVE :: aux3 !<
109	REAL(wp), DIMENSION(nau1), SAVE :: auy3 !<
110
111	#elif defined( __nec_fft )
112	INTEGER(iwp), SAVE :: nz1 !<
113
114	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: trig_xb !<
115	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: trig_xf !<
116	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: trig_yb !<
117	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: trig_yf !<
118
119	#elif defined( __cuda_fft )
120	INTEGER(C_INT), SAVE :: plan_xf !<
121	INTEGER(C_INT), SAVE :: plan_xi !<
122	INTEGER(C_INT), SAVE :: plan_yf !<
123	INTEGER(C_INT), SAVE :: plan_yi !<
124
125	#endif
126
127	#if defined( __fftw )
128	INCLUDE 'fftw3.f03'
129	COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE :: x_out !<
130	COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE :: y_out !<
131
132	INTEGER(KIND=C_INT) :: nx_c !<
133	INTEGER(KIND=C_INT) :: ny_c !<
134
135	REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE :: x_in !<
136	REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE :: y_in !<
137	!$OMP THREADPRIVATE( x_out, y_out, x_in, y_in )
138
139
140	TYPE(C_PTR), SAVE :: plan_xf, plan_xi, plan_yf, plan_yi
141	#endif
142
143	!
144	!-- Public interfaces
145	INTERFACE fft_init
146	MODULE PROCEDURE fft_init
147	END INTERFACE fft_init
148
149	INTERFACE fft_x
150	MODULE PROCEDURE fft_x
151	END INTERFACE fft_x
152
153	INTERFACE fft_x_1d
154	MODULE PROCEDURE fft_x_1d
155	END INTERFACE fft_x_1d
156
157	INTERFACE fft_y
158	MODULE PROCEDURE fft_y
159	END INTERFACE fft_y
160
161	INTERFACE fft_y_1d
162	MODULE PROCEDURE fft_y_1d
163	END INTERFACE fft_y_1d
164
165	INTERFACE fft_x_m
166	MODULE PROCEDURE fft_x_m
167	END INTERFACE fft_x_m
168
169	INTERFACE fft_y_m
170	MODULE PROCEDURE fft_y_m
171	END INTERFACE fft_y_m
172
173	CONTAINS
174
175
176	!--------------------------------------------------------------------------------------------------!
177	! Description:
178	! ------------
179	!> @todo Missing subroutine description.
180	!--------------------------------------------------------------------------------------------------!
181	SUBROUTINE fft_init
182
183	USE pegrid, &
184	ONLY: pdims
185
186	IMPLICIT NONE
187
188	!
189	!-- The following temporary working arrays have to be on stack or private
190	!-- in OpenMP sense
191	#if defined( __ibm )
192	REAL(wp), DIMENSION(nau2) :: aux2 !<
193	REAL(wp), DIMENSION(nau2) :: auy2 !<
194	REAL(wp), DIMENSION(nau2) :: aux4 !<
195	REAL(wp), DIMENSION(nau2) :: auy4 !<
196	REAL(wp), DIMENSION(0:nx+2) :: workx !<
197	REAL(wp), DIMENSION(0:ny+2) :: worky !<
198	#elif defined( __nec_fft )
199	REAL(wp), DIMENSION(0:nx+3,nz+1) :: work_x !<
200	REAL(wp), DIMENSION(0:ny+3,nz+1) :: work_y !<
201	REAL(wp), DIMENSION(6*(nx+3),nz+1) :: workx !<
202	REAL(wp), DIMENSION(6*(ny+3),nz+1) :: worky !<
203	#endif
204
205	!
206	!-- Return, if already called
207	IF ( init_fft ) THEN
208	RETURN
209	ELSE
210	init_fft = .TRUE.
211	ENDIF
212
213	#if defined( _OPENACC ) && defined( __cuda_fft )
214	fft_method = 'system-specific'
215	#endif
216
217	!
218	!-- Switch to tell the Poisson-solver that the vectorized version of Temperton-fft is to be used.
219	IF ( fft_method == 'temperton-algorithm' .AND. loop_optimization == 'vector' .AND. &
220	pdims(1) /= 1 .AND. pdims(2) /= 1 ) THEN
221	temperton_fft_vec = .TRUE.
222	ENDIF
223
224	IF ( fft_method == 'system-specific' ) THEN
225
226	dnx = 1.0_wp / ( nx + 1.0_wp )
227	dny = 1.0_wp / ( ny + 1.0_wp )
228	sqr_dnx = SQRT( dnx )
229	sqr_dny = SQRT( dny )
230	#if defined( __ibm )
231	!
232	!-- Initialize tables for fft along x
233	CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, aux2, nau2 )
234	CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, aux4, nau2 )
235	!
236	!-- Initialize tables for fft along y
237	CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, auy2, nau2 )
238	CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, auy4, nau2 )
239	#elif defined( __nec_fft )
240	message_string = 'fft method "' // TRIM( fft_method) // '" currently does not work on NEC'
241	CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 )
242
243	ALLOCATE( trig_xb(2(nx+1)), trig_xf(2(nx+1)), trig_yb(2(ny+1)), trig_yf(2(ny+1)) )
244
245	work_x = 0.0_wp
246	work_y = 0.0_wp
247	nz1 = nz + MOD( nz+1, 2 ) ! odd nz slows down fft significantly
248	! when using the NEC ffts
249
250	!
251	!-- Initialize tables for fft along x (non-vector and vector case (M))
252	CALL DZFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xf, workx, 0 )
253	CALL ZDFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xb, workx, 0 )
254	CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, trig_xf, workx, 0 )
255	CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, trig_xb, workx, 0 )
256	!
257	!-- Initialize tables for fft along y (non-vector and vector case (M))
258	CALL DZFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yf, worky, 0 )
259	CALL ZDFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yb, worky, 0 )
260	CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, trig_yf, worky, 0 )
261	CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, trig_yb, worky, 0 )
262	#elif defined( __cuda_fft )
263	CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
264	CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
265	CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
266	CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
267	#else
268	message_string = 'no system-specific fft-call available'
269	CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
270	#endif
271	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
272	!
273	!-- Temperton-algorithm
274	!-- Initialize tables for fft along x and y
275	ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
276
277	CALL set99( trigs_x, ifax_x, nx+1 )
278	CALL set99( trigs_y, ifax_y, ny+1 )
279
280	IF ( temperton_fft_vec ) THEN
281	ALLOCATE( f_vec_x((nyn_x-nys_x+1)*(nzt_x-nzb_x+1),0:nx+2) )
282	ENDIF
283
284
285
286	ELSEIF ( fft_method == 'fftw' ) THEN
287	!
288	!-- FFTW
289	#if defined( __fftw )
290	nx_c = nx+1
291	ny_c = ny+1
292	!$OMP PARALLEL
293	ALLOCATE( x_in(0:nx+2), y_in(0:ny+2), x_out(0:(nx+1)/2), y_out(0:(ny+1)/2) )
294	!$OMP END PARALLEL
295	plan_xf = FFTW_PLAN_DFT_R2C_1D( nx_c, x_in, x_out, FFTW_ESTIMATE )
296	plan_xi = FFTW_PLAN_DFT_C2R_1D( nx_c, x_out, x_in, FFTW_ESTIMATE )
297	plan_yf = FFTW_PLAN_DFT_R2C_1D( ny_c, y_in, y_out, FFTW_ESTIMATE )
298	plan_yi = FFTW_PLAN_DFT_C2R_1D( ny_c, y_out, y_in, FFTW_ESTIMATE )
299	#else
300	message_string = 'preprocessor switch for fftw is missing'
301	CALL message( 'fft_init', 'PA0080', 1, 2, 0, 6, 0 )
302	#endif
303
304	ELSEIF ( fft_method == 'singleton-algorithm' ) THEN
305
306	CONTINUE
307
308	ELSE
309
310	message_string = 'fft method "' // TRIM( fft_method) // '" not available'
311	CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
312	ENDIF
313
314	END SUBROUTINE fft_init
315
316
317	!--------------------------------------------------------------------------------------------------!
318	! Description:
319	! ------------
320	!> Fourier-transformation along x-direction.
321	!> Version for 2D-decomposition.
322	!> It uses internal algorithms (Singleton or Temperton) or system-specific routines, if they are
323	!> available.
324	!--------------------------------------------------------------------------------------------------!
325
326	SUBROUTINE fft_x( ar, direction, ar_2d, ar_inv )
327
328
329	IMPLICIT NONE
330
331	CHARACTER (LEN=*) :: direction !<
332
333	COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: cwork !<
334
335	INTEGER(iwp) :: i !<
336	INTEGER(iwp) :: ishape(1) !<
337	INTEGER(iwp) :: j !<
338	INTEGER(iwp) :: k !<
339	INTEGER(iwp) :: mm !<
340
341	LOGICAL :: forward_fft !<
342
343	REAL(wp), DIMENSION(0:nx+2) :: work !<
344	REAL(wp), DIMENSION(nx+2) :: work1 !<
345
346	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: work_vec !<
347	REAL(wp), DIMENSION(0:nx,nys_x:nyn_x), OPTIONAL :: ar_2d !<
348
349	REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: ar !<
350	REAL(wp), DIMENSION(nys_x:nyn_x,nzb_x:nzt_x,0:nx), OPTIONAL :: ar_inv !<
351
352	#if defined( __ibm )
353	REAL(wp), DIMENSION(nau2) :: aux2 !<
354	REAL(wp), DIMENSION(nau2) :: aux4 !<
355	#elif defined( __nec_fft )
356	REAL(wp), DIMENSION(6*(nx+1)) :: work2 !<
357	#elif defined( __cuda_fft )
358	COMPLEX(dp), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) :: ar_tmp !<
359	!$ACC DECLARE CREATE(ar_tmp)
360	#endif
361
362	!
363	!-- To avoid compiler warning: Unused dummy argument âar_2dâ
364	IF ( PRESENT( ar_2d ) ) CONTINUE
365
366	IF ( direction == 'forward' ) THEN
367	forward_fft = .TRUE.
368	ELSE
369	forward_fft = .FALSE.
370	ENDIF
371
372	IF ( fft_method == 'singleton-algorithm' ) THEN
373
374	!
375	!-- Performing the fft with singleton's software works on every system, since it is part of
376	!-- the model.
377	ALLOCATE( cwork(0:nx) )
378
379	IF ( forward_fft ) THEN
380
381	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
382	!$OMP DO
383	DO k = nzb_x, nzt_x
384	DO j = nys_x, nyn_x
385
386	DO i = 0, nx
387	cwork(i) = CMPLX( ar(i,j,k), KIND=wp )
388	ENDDO
389
390	ishape = SHAPE( cwork )
391	CALL FFTN( cwork, ishape )
392
393	DO i = 0, (nx+1)/2
394	ar(i,j,k) = REAL( cwork(i), KIND=wp )
395	ENDDO
396	DO i = 1, (nx+1)/2 - 1
397	ar(nx+1-i,j,k) = -AIMAG( cwork(i) )
398	ENDDO
399
400	ENDDO
401	ENDDO
402	!$OMP END PARALLEL
403
404	ELSE
405
406	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
407	!$OMP DO
408	DO k = nzb_x, nzt_x
409	DO j = nys_x, nyn_x
410
411	cwork(0) = CMPLX( ar(0,j,k), 0.0_wp, KIND=wp )
412	DO i = 1, (nx+1)/2 - 1
413	cwork(i) = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k), KIND=wp )
414	cwork(nx+1-i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k), KIND=wp )
415	ENDDO
416	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp, KIND=wp )
417
418	ishape = SHAPE( cwork )
419	CALL FFTN( cwork, ishape, inv = .TRUE. )
420
421	DO i = 0, nx
422	ar(i,j,k) = REAL( cwork(i), KIND=wp )
423	ENDDO
424
425	ENDDO
426	ENDDO
427	!$OMP END PARALLEL
428
429	ENDIF
430
431	DEALLOCATE( cwork )
432
433	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
434
435	!
436	!-- Performing the fft with Temperton's software works on every system, since it is part of
437	!-- the model.
438	IF ( forward_fft ) THEN
439
440	IF ( .NOT. temperton_fft_vec ) THEN
441
442	!$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
443	!$OMP DO
444	DO k = nzb_x, nzt_x
445	DO j = nys_x, nyn_x
446
447	work(0:nx) = ar(0:nx,j,k)
448	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
449
450	DO i = 0, (nx+1)/2
451	ar(i,j,k) = work(2*i)
452	ENDDO
453	DO i = 1, (nx+1)/2 - 1
454	ar(nx+1-i,j,k) = work(2*i+1)
455	ENDDO
456
457	ENDDO
458	ENDDO
459	!$OMP END PARALLEL
460
461	ELSE
462
463	!
464	!-- Vector version of the Temperton-algorithm. Computes multiple 1-D FFT's.
465	ALLOCATE( work_vec( (nyn_x-nys_x+1)*(nzt_x-nzb_x+1),nx+2) )
466	!
467	!-- f_vec_x is already set in transpose_zx
468	CALL fft991cy_vec( f_vec_x, work_vec, trigs_x, ifax_x, nx+1, -1 )
469	DEALLOCATE( work_vec )
470
471	IF ( PRESENT( ar_inv ) ) THEN
472
473	DO k = nzb_x, nzt_x
474	DO j = nys_x, nyn_x
475	mm = j-nys_x+1+(k-nzb_x)*(nyn_x-nys_x+1)
476	DO i = 0, (nx+1)/2
477	ar_inv(j,k,i) = f_vec_x(mm,2*i)
478	ENDDO
479	DO i = 1, (nx+1)/2-1
480	ar_inv(j,k,nx+1-i) = f_vec_x(mm,2*i+1)
481	ENDDO
482	ENDDO
483	ENDDO
484
485	ELSE
486
487	DO k = nzb_x, nzt_x
488	DO j = nys_x, nyn_x
489	mm = j-nys_x+1+(k-nzb_x)*(nyn_x-nys_x+1)
490	DO i = 0, (nx+1)/2
491	ar(i,j,k) = f_vec_x(mm,2*i)
492	ENDDO
493	DO i = 1, (nx+1)/2-1
494	ar(nx+1-i,j,k) = f_vec_x(mm,2*i+1)
495	ENDDO
496	ENDDO
497	ENDDO
498
499	ENDIF
500
501	ENDIF
502
503	ELSE
504
505	!
506	!-- Backward fft
507	IF ( .NOT. temperton_fft_vec ) THEN
508
509	!$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
510	!$OMP DO
511	DO k = nzb_x, nzt_x
512	DO j = nys_x, nyn_x
513
514	DO i = 0, (nx+1)/2
515	work(2*i) = ar(i,j,k)
516	ENDDO
517	DO i = 1, (nx+1)/2 - 1
518	work(2*i+1) = ar(nx+1-i,j,k)
519	ENDDO
520	work(1) = 0.0_wp
521	work(nx+2) = 0.0_wp
522
523	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
524	ar(0:nx,j,k) = work(0:nx)
525
526	ENDDO
527	ENDDO
528	!$OMP END PARALLEL
529
530	ELSE
531
532	IF ( PRESENT( ar_inv ) ) THEN
533
534	DO k = nzb_x, nzt_x
535	DO j = nys_x, nyn_x
536	mm = j-nys_x+1+(k-nzb_x)*(nyn_x-nys_x+1)
537	DO i = 0, (nx+1)/2
538	f_vec_x(mm,2*i) = ar_inv(j,k,i)
539	ENDDO
540	DO i = 1, (nx+1)/2-1
541	f_vec_x(mm,2*i+1) = ar_inv(j,k,nx+1-i)
542	ENDDO
543	ENDDO
544	ENDDO
545
546	ELSE
547
548	DO k = nzb_x, nzt_x
549	DO j = nys_x, nyn_x
550	mm = j-nys_x+1+(k-nzb_x)*(nyn_x-nys_x+1)
551	DO i = 0, (nx+1)/2
552	f_vec_x(mm,2*i) = ar(i,j,k)
553	ENDDO
554	DO i = 1, (nx+1)/2-1
555	f_vec_x(mm,2*i+1) = ar(nx+1-i,j,k)
556	ENDDO
557	ENDDO
558	ENDDO
559
560	ENDIF
561	f_vec_x(:,1) = 0.0_wp
562	f_vec_x(:,nx+2) = 0.0_wp
563
564	ALLOCATE( work_vec((nyn_x-nys_x+1)*(nzt_x-nzb_x+1),nx+2) )
565	CALL fft991cy_vec( f_vec_x, work_vec, trigs_x, ifax_x, nx+1, 1 )
566	DEALLOCATE( work_vec )
567
568	ENDIF
569
570	ENDIF
571
572	ELSEIF ( fft_method == 'fftw' ) THEN
573
574	#if defined( __fftw )
575	IF ( forward_fft ) THEN
576
577	!$OMP PARALLEL PRIVATE ( work, i, j, k )
578	!$OMP DO
579	DO k = nzb_x, nzt_x
580	DO j = nys_x, nyn_x
581
582	x_in(0:nx) = ar(0:nx,j,k)
583	CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
584
585	IF ( PRESENT( ar_2d ) ) THEN
586
587	DO i = 0, (nx+1)/2
588	ar_2d(i,j) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
589	ENDDO
590	DO i = 1, (nx+1)/2 - 1
591	ar_2d(nx+1-i,j) = AIMAG( x_out(i) ) / ( nx+1 )
592	ENDDO
593
594	ELSE
595
596	DO i = 0, (nx+1)/2
597	ar(i,j,k) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
598	ENDDO
599	DO i = 1, (nx+1)/2 - 1
600	ar(nx+1-i,j,k) = AIMAG( x_out(i) ) / ( nx+1 )
601	ENDDO
602
603	ENDIF
604
605	ENDDO
606	ENDDO
607	!$OMP END PARALLEL
608
609	ELSE
610	!$OMP PARALLEL PRIVATE ( work, i, j, k )
611	!$OMP DO
612	DO k = nzb_x, nzt_x
613	DO j = nys_x, nyn_x
614
615	IF ( PRESENT( ar_2d ) ) THEN
616
617	x_out(0) = CMPLX( ar_2d(0,j), 0.0_wp, KIND=wp )
618	DO i = 1, (nx+1)/2 - 1
619	x_out(i) = CMPLX( ar_2d(i,j), ar_2d(nx+1-i,j), KIND=wp )
620	ENDDO
621	x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0_wp, KIND=wp )
622
623	ELSE
624
625	x_out(0) = CMPLX( ar(0,j,k), 0.0_wp, KIND=wp )
626	DO i = 1, (nx+1)/2 - 1
627	x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k), KIND=wp )
628	ENDDO
629	x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp, KIND=wp )
630
631	ENDIF
632
633	CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
634	ar(0:nx,j,k) = x_in(0:nx)
635
636	ENDDO
637	ENDDO
638	!$OMP END PARALLEL
639
640	ENDIF
641	#endif
642
643	ELSEIF ( fft_method == 'system-specific' ) THEN
644
645	#if defined( __ibm )
646	IF ( forward_fft ) THEN
647
648	!$OMP PARALLEL PRIVATE ( work, i, j, k )
649	!$OMP DO
650	DO k = nzb_x, nzt_x
651	DO j = nys_x, nyn_x
652
653	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, aux2, nau2 )
654
655	DO i = 0, (nx+1)/2
656	ar(i,j,k) = work(2*i)
657	ENDDO
658	DO i = 1, (nx+1)/2 - 1
659	ar(nx+1-i,j,k) = work(2*i+1)
660	ENDDO
661
662	ENDDO
663	ENDDO
664	!$OMP END PARALLEL
665
666	ELSE
667
668	!$OMP PARALLEL PRIVATE ( work, i, j, k )
669	!$OMP DO
670	DO k = nzb_x, nzt_x
671	DO j = nys_x, nyn_x
672
673	DO i = 0, (nx+1)/2
674	work(2*i) = ar(i,j,k)
675	ENDDO
676	DO i = 1, (nx+1)/2 - 1
677	work(2*i+1) = ar(nx+1-i,j,k)
678	ENDDO
679	work(1) = 0.0_wp
680	work(nx+2) = 0.0_wp
681
682	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, aux4, nau2 )
683
684	DO i = 0, nx
685	ar(i,j,k) = work(i)
686	ENDDO
687
688	ENDDO
689	ENDDO
690	!$OMP END PARALLEL
691
692	ENDIF
693
694	#elif defined( __nec_fft )
695
696	IF ( forward_fft ) THEN
697
698	!$OMP PARALLEL PRIVATE ( work, i, j, k )
699	!$OMP DO
700	DO k = nzb_x, nzt_x
701	DO j = nys_x, nyn_x
702
703	work(0:nx) = ar(0:nx,j,k)
704
705	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
706
707	DO i = 0, (nx+1)/2
708	ar(i,j,k) = work(2*i)
709	ENDDO
710	DO i = 1, (nx+1)/2 - 1
711	ar(nx+1-i,j,k) = work(2*i+1)
712	ENDDO
713
714	ENDDO
715	ENDDO
716	!$END OMP PARALLEL
717
718	ELSE
719
720	!$OMP PARALLEL PRIVATE ( work, i, j, k )
721	!$OMP DO
722	DO k = nzb_x, nzt_x
723	DO j = nys_x, nyn_x
724
725	DO i = 0, (nx+1)/2
726	work(2*i) = ar(i,j,k)
727	ENDDO
728	DO i = 1, (nx+1)/2 - 1
729	work(2*i+1) = ar(nx+1-i,j,k)
730	ENDDO
731	work(1) = 0.0_wp
732	work(nx+2) = 0.0_wp
733
734	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
735
736	ar(0:nx,j,k) = work(0:nx)
737
738	ENDDO
739	ENDDO
740	!$OMP END PARALLEL
741
742	ENDIF
743
744	#elif defined( __cuda_fft )
745
746	IF ( forward_fft ) THEN
747
748	!$ACC HOST_DATA USE_DEVICE(ar, ar_tmp)
749	CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
750	!$ACC END HOST_DATA
751
752	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i,j,k) &
753	!$ACC PRESENT(ar, ar_tmp)
754	DO k = nzb_x, nzt_x
755	DO j = nys_x, nyn_x
756
757	DO i = 0, (nx+1)/2
758	ar(i,j,k) = REAL( ar_tmp(i,j,k), KIND=wp ) * dnx
759	ENDDO
760
761	DO i = 1, (nx+1)/2 - 1
762	ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
763	ENDDO
764
765	ENDDO
766	ENDDO
767
768	ELSE
769
770	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i,j,k) &
771	!$ACC PRESENT(ar, ar_tmp)
772	DO k = nzb_x, nzt_x
773	DO j = nys_x, nyn_x
774
775	ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0_wp, KIND=wp )
776
777	DO i = 1, (nx+1)/2 - 1
778	ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k), KIND=wp )
779	ENDDO
780	ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp, KIND=wp )
781
782	ENDDO
783	ENDDO
784
785	!$ACC HOST_DATA USE_DEVICE(ar, ar_tmp)
786	CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
787	!$ACC END HOST_DATA
788
789	ENDIF
790
791	#endif
792
793	ENDIF
794
795	END SUBROUTINE fft_x
796
797	!--------------------------------------------------------------------------------------------------!
798	! Description:
799	! ------------
800	!> Fourier-transformation along x-direction.
801	!> Version for 1D-decomposition.
802	!> It uses internal algorithms (Singleton or Temperton) or system-specific routines, if they are
803	!> available.
804	!--------------------------------------------------------------------------------------------------!
805
806	SUBROUTINE fft_x_1d( ar, direction )
807
808
809	IMPLICIT NONE
810
811	CHARACTER (LEN=*) :: direction !<
812
813	INTEGER(iwp) :: i !<
814	INTEGER(iwp) :: ishape(1) !<
815
816	LOGICAL :: forward_fft !<
817
818	REAL(wp), DIMENSION(0:nx) :: ar !<
819	REAL(wp), DIMENSION(0:nx+2) :: work !<
820	REAL(wp), DIMENSION(nx+2) :: work1 !<
821
822	COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: cwork !<
823
824	#if defined( __ibm )
825	REAL(wp), DIMENSION(nau2) :: aux2 !<
826	REAL(wp), DIMENSION(nau2) :: aux4 !<
827	#elif defined( __nec_fft )
828	REAL(wp), DIMENSION(6*(nx+1)) :: work2 !<
829	#endif
830
831	IF ( direction == 'forward' ) THEN
832	forward_fft = .TRUE.
833	ELSE
834	forward_fft = .FALSE.
835	ENDIF
836
837	IF ( fft_method == 'singleton-algorithm' ) THEN
838
839	!
840	!-- Performing the fft with singleton's software works on every system, since it is part of
841	!-- the model.
842	ALLOCATE( cwork(0:nx) )
843
844	IF ( forward_fft ) THEN
845
846	DO i = 0, nx
847	cwork(i) = CMPLX( ar(i), KIND=wp )
848	ENDDO
849	ishape = SHAPE( cwork )
850	CALL FFTN( cwork, ishape )
851	DO i = 0, (nx+1)/2
852	ar(i) = REAL( cwork(i), KIND=wp )
853	ENDDO
854	DO i = 1, (nx+1)/2 - 1
855	ar(nx+1-i) = -AIMAG( cwork(i) )
856	ENDDO
857
858	ELSE
859
860	cwork(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
861	DO i = 1, (nx+1)/2 - 1
862	cwork(i) = CMPLX( ar(i), -ar(nx+1-i), KIND=wp )
863	cwork(nx+1-i) = CMPLX( ar(i), ar(nx+1-i), KIND=wp )
864	ENDDO
865	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp, KIND=wp )
866
867	ishape = SHAPE( cwork )
868	CALL FFTN( cwork, ishape, inv = .TRUE. )
869
870	DO i = 0, nx
871	ar(i) = REAL( cwork(i), KIND=wp )
872	ENDDO
873
874	ENDIF
875
876	DEALLOCATE( cwork )
877
878	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
879
880	!
881	!-- Performing the fft with Temperton's software works on every system, since it is part of
882	!-- the model.
883	IF ( forward_fft ) THEN
884
885	work(0:nx) = ar
886	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
887
888	DO i = 0, (nx+1)/2
889	ar(i) = work(2*i)
890	ENDDO
891	DO i = 1, (nx+1)/2 - 1
892	ar(nx+1-i) = work(2*i+1)
893	ENDDO
894
895	ELSE
896
897	DO i = 0, (nx+1)/2
898	work(2*i) = ar(i)
899	ENDDO
900	DO i = 1, (nx+1)/2 - 1
901	work(2*i+1) = ar(nx+1-i)
902	ENDDO
903	work(1) = 0.0_wp
904	work(nx+2) = 0.0_wp
905
906	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
907	ar = work(0:nx)
908
909	ENDIF
910
911	ELSEIF ( fft_method == 'fftw' ) THEN
912
913	#if defined( __fftw )
914	IF ( forward_fft ) THEN
915
916	x_in(0:nx) = ar(0:nx)
917	CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
918
919	DO i = 0, (nx+1)/2
920	ar(i) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
921	ENDDO
922	DO i = 1, (nx+1)/2 - 1
923	ar(nx+1-i) = AIMAG( x_out(i) ) / ( nx+1 )
924	ENDDO
925
926	ELSE
927
928	x_out(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
929	DO i = 1, (nx+1)/2 - 1
930	x_out(i) = CMPLX( ar(i), ar(nx+1-i), KIND=wp )
931	ENDDO
932	x_out((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp, KIND=wp )
933
934	CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
935	ar(0:nx) = x_in(0:nx)
936
937	ENDIF
938	#endif
939
940	ELSEIF ( fft_method == 'system-specific' ) THEN
941
942	#if defined( __ibm )
943	IF ( forward_fft ) THEN
944
945	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, aux2, nau2 )
946
947	DO i = 0, (nx+1)/2
948	ar(i) = work(2*i)
949	ENDDO
950	DO i = 1, (nx+1)/2 - 1
951	ar(nx+1-i) = work(2*i+1)
952	ENDDO
953
954	ELSE
955
956	DO i = 0, (nx+1)/2
957	work(2*i) = ar(i)
958	ENDDO
959	DO i = 1, (nx+1)/2 - 1
960	work(2*i+1) = ar(nx+1-i)
961	ENDDO
962	work(1) = 0.0_wp
963	work(nx+2) = 0.0_wp
964
965	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, aux4, nau2 )
966
967	DO i = 0, nx
968	ar(i) = work(i)
969	ENDDO
970
971	ENDIF
972	#elif defined( __nec_fft )
973	IF ( forward_fft ) THEN
974
975	work(0:nx) = ar(0:nx)
976
977	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
978
979	DO i = 0, (nx+1)/2
980	ar(i) = work(2*i)
981	ENDDO
982	DO i = 1, (nx+1)/2 - 1
983	ar(nx+1-i) = work(2*i+1)
984	ENDDO
985
986	ELSE
987
988	DO i = 0, (nx+1)/2
989	work(2*i) = ar(i)
990	ENDDO
991	DO i = 1, (nx+1)/2 - 1
992	work(2*i+1) = ar(nx+1-i)
993	ENDDO
994	work(1) = 0.0_wp
995	work(nx+2) = 0.0_wp
996
997	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
998
999	ar(0:nx) = work(0:nx)
1000
1001	ENDIF
1002	#endif
1003
1004	ENDIF
1005
1006	END SUBROUTINE fft_x_1d
1007
1008	!--------------------------------------------------------------------------------------------------!
1009	! Description:
1010	! ------------
1011	!> Fourier-transformation along y-direction.
1012	!> Version for 2D-decomposition.
1013	!> It uses internal algorithms (Singleton or Temperton) or system-specific routines, if they are
1014	!> available.
1015	!>
1016	!> direction: 'forward' or 'backward'
1017	!> ar, ar_tr: 3D data arrays
1018	!> forward: ar: before ar_tr: after transformation
1019	!> backward: ar_tr: before ar: after transfosition
1020	!>
1021	!> In case of non-overlapping transposition/transformation:
1022	!> nxl_y_bound = nxl_y_l = nxl_y
1023	!> nxr_y_bound = nxr_y_l = nxr_y
1024	!>
1025	!> In case of overlapping transposition/transformation
1026	!> - nxl_y_bound and nxr_y_bound have the original values of nxl_y, nxr_y. ar_tr is dimensioned
1027	!> using these values.
1028	!> - nxl_y_l = nxr_y_r. ar is dimensioned with these values, so that transformation is carried out
1029	!> for a 2D-plane only.
1030	!--------------------------------------------------------------------------------------------------!
1031
1032	SUBROUTINE fft_y( ar, direction, ar_tr, nxl_y_bound, nxr_y_bound, nxl_y_l, nxr_y_l, ar_inv )
1033
1034
1035	IMPLICIT NONE
1036
1037	CHARACTER (LEN=*) :: direction !<
1038
1039	INTEGER(iwp) :: i !<
1040	INTEGER(iwp) :: j !<
1041	INTEGER(iwp) :: jshape(1) !<
1042	INTEGER(iwp) :: k !<
1043	INTEGER(iwp) :: mm !<
1044	INTEGER(iwp) :: nxl_y_bound !<
1045	INTEGER(iwp) :: nxl_y_l !<
1046	INTEGER(iwp) :: nxr_y_bound !<
1047	INTEGER(iwp) :: nxr_y_l !<
1048
1049	LOGICAL :: forward_fft !<
1050
1051	REAL(wp), DIMENSION(0:ny+2) :: work !<
1052	REAL(wp), DIMENSION(ny+2) :: work1 !<
1053
1054	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: f_vec_y
1055	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: work_vec
1056
1057	REAL(wp), DIMENSION(0:ny,nxl_y_l:nxr_y_l,nzb_y:nzt_y) :: ar !<
1058	REAL(wp), DIMENSION(nxl_y:nxr_y,nzb_y:nzt_y,0:ny), OPTIONAL :: ar_inv !<
1059	REAL(wp), DIMENSION(0:ny,nxl_y_bound:nxr_y_bound,nzb_y:nzt_y), OPTIONAL :: ar_tr !<
1060
1061	COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: cwork !<
1062
1063	#if defined( __ibm )
1064	REAL(wp), DIMENSION(nau2) :: auy2 !<
1065	REAL(wp), DIMENSION(nau2) :: auy4 !<
1066	#elif defined( __nec_fft )
1067	REAL(wp), DIMENSION(6*(ny+1)) :: work2 !<
1068	#elif defined( __cuda_fft )
1069	COMPLEX(dp), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) :: ar_tmp !<
1070	!$ACC DECLARE CREATE(ar_tmp)
1071	#endif
1072
1073
1074	IF ( direction == 'forward' ) THEN
1075	forward_fft = .TRUE.
1076	ELSE
1077	forward_fft = .FALSE.
1078	ENDIF
1079
1080	IF ( fft_method == 'singleton-algorithm' ) THEN
1081
1082	!
1083	!-- Performing the fft with singleton's software works on every system, since it is part of
1084	!-- the model.
1085	ALLOCATE( cwork(0:ny) )
1086
1087	IF ( forward_fft ) THEN
1088
1089	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
1090	!$OMP DO
1091	DO k = nzb_y, nzt_y
1092	DO i = nxl_y_l, nxr_y_l
1093
1094	DO j = 0, ny
1095	cwork(j) = CMPLX( ar(j,i,k), KIND=wp )
1096	ENDDO
1097
1098	jshape = SHAPE( cwork )
1099	CALL FFTN( cwork, jshape )
1100
1101	DO j = 0, (ny+1)/2
1102	ar_tr(j,i,k) = REAL( cwork(j), KIND=wp )
1103	ENDDO
1104	DO j = 1, (ny+1)/2 - 1
1105	ar_tr(ny+1-j,i,k) = -AIMAG( cwork(j) )
1106	ENDDO
1107
1108	ENDDO
1109	ENDDO
1110	!$OMP END PARALLEL
1111
1112	ELSE
1113
1114	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
1115	!$OMP DO
1116	DO k = nzb_y, nzt_y
1117	DO i = nxl_y_l, nxr_y_l
1118
1119	cwork(0) = CMPLX( ar_tr(0,i,k), 0.0_wp, KIND=wp )
1120	DO j = 1, (ny+1)/2 - 1
1121	cwork(j) = CMPLX( ar_tr(j,i,k), -ar_tr(ny+1-j,i,k), KIND=wp )
1122	cwork(ny+1-j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k), KIND=wp )
1123	ENDDO
1124	cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp, KIND=wp )
1125
1126	jshape = SHAPE( cwork )
1127	CALL FFTN( cwork, jshape, inv = .TRUE. )
1128
1129	DO j = 0, ny
1130	ar(j,i,k) = REAL( cwork(j), KIND=wp )
1131	ENDDO
1132
1133	ENDDO
1134	ENDDO
1135	!$OMP END PARALLEL
1136
1137	ENDIF
1138
1139	DEALLOCATE( cwork )
1140
1141	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
1142
1143	!
1144	!-- Performing the fft with Temperton's software works on every system, since it is part of
1145	!-- the model.
1146	IF ( forward_fft ) THEN
1147
1148	IF ( .NOT. temperton_fft_vec ) THEN
1149
1150	!$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
1151	!$OMP DO
1152	DO k = nzb_y, nzt_y
1153	DO i = nxl_y_l, nxr_y_l
1154
1155	work(0:ny) = ar(0:ny,i,k)
1156	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
1157
1158	DO j = 0, (ny+1)/2
1159	ar_tr(j,i,k) = work(2*j)
1160	ENDDO
1161	DO j = 1, (ny+1)/2 - 1
1162	ar_tr(ny+1-j,i,k) = work(2*j+1)
1163	ENDDO
1164
1165	ENDDO
1166	ENDDO
1167	!$OMP END PARALLEL
1168
1169	ELSE
1170	!
1171	!-- Vector version of Temperton-fft. Computes multiple 1-D FFT's.
1172	ALLOCATE( f_vec_y((nxr_y_l-nxl_y_l+1)*(nzt_y-nzb_y+1),0:ny+2) )
1173
1174	mm = 1
1175	DO k = nzb_y, nzt_y
1176	DO i = nxl_y_l, nxr_y_l
1177	f_vec_y(mm,0:nx) = ar(0:nx,i,k)
1178	mm = mm+1
1179	ENDDO
1180	ENDDO
1181
1182	ALLOCATE( work_vec( (nxr_y_l-nxl_y_l+1)*(nzt_y-nzb_y+1),ny+2) )
1183	CALL fft991cy_vec( f_vec_y, work_vec, trigs_y, ifax_y, ny+1, -1 )
1184	DEALLOCATE( work_vec )
1185
1186	IF( PRESENT( ar_inv ) ) THEN
1187
1188	DO k = nzb_y, nzt_y
1189	DO i = nxl_y_l, nxr_y_l
1190	mm = i-nxl_y_l+1+(k-nzb_y)*(nxr_y_l-nxl_y_l+1)
1191	DO j = 0, (ny+1)/2
1192	ar_inv(i,k,j) = f_vec_y(mm,2*j)
1193	ENDDO
1194	DO j = 1, (ny+1)/2 - 1
1195	ar_inv(i,k,ny+1-j) = f_vec_y(mm,2*j+1)
1196	ENDDO
1197	ENDDO
1198	ENDDO
1199
1200	ELSE
1201
1202	DO k = nzb_y, nzt_y
1203	DO i = nxl_y_l, nxr_y_l
1204	mm = i-nxl_y_l+1+(k-nzb_y)*(nxr_y_l-nxl_y_l+1)
1205	DO j = 0, (ny+1)/2
1206	ar(j,i,k) = f_vec_y(mm,2*j)
1207	ENDDO
1208	DO j = 1, (ny+1)/2 - 1
1209	ar(ny+1-j,i,k) = f_vec_y(mm,2*j+1)
1210	ENDDO
1211	ENDDO
1212	ENDDO
1213
1214	ENDIF
1215
1216	DEALLOCATE( f_vec_y )
1217
1218	ENDIF
1219
1220	ELSE
1221
1222	IF ( .NOT. temperton_fft_vec ) THEN
1223
1224	!$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
1225	!$OMP DO
1226	DO k = nzb_y, nzt_y
1227	DO i = nxl_y_l, nxr_y_l
1228
1229	DO j = 0, (ny+1)/2
1230	work(2*j) = ar_tr(j,i,k)
1231	ENDDO
1232	DO j = 1, (ny+1)/2 - 1
1233	work(2*j+1) = ar_tr(ny+1-j,i,k)
1234	ENDDO
1235	work(1) = 0.0_wp
1236	work(ny+2) = 0.0_wp
1237
1238	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
1239	ar(0:ny,i,k) = work(0:ny)
1240
1241	ENDDO
1242	ENDDO
1243	!$OMP END PARALLEL
1244
1245	ELSE
1246
1247	ALLOCATE( f_vec_y((nxr_y_l-nxl_y_l+1)*(nzt_y-nzb_y+1),0:ny+2) )
1248
1249	IF ( PRESENT( ar_inv ) ) THEN
1250
1251	DO k = nzb_y, nzt_y
1252	DO i = nxl_y_l, nxr_y_l
1253	mm = i-nxl_y_l+1+(k-nzb_y)*(nxr_y_l-nxl_y_l+1)
1254	DO j = 0, (ny+1)/2
1255	f_vec_y(mm,2*j) = ar_inv(i,k,j)
1256	ENDDO
1257	DO j = 1, (ny+1)/2 - 1
1258	f_vec_y(mm,2*j+1) = ar_inv(i,k,ny+1-j)
1259	ENDDO
1260	ENDDO
1261	ENDDO
1262
1263	ELSE
1264
1265	DO k = nzb_y, nzt_y
1266	DO i = nxl_y_l, nxr_y_l
1267	mm = i-nxl_y_l+1+(k-nzb_y)*(nxr_y_l-nxl_y_l+1)
1268	DO j = 0, (ny+1)/2
1269	f_vec_y(mm,2*j) = ar(j,i,k)
1270	ENDDO
1271	DO j = 1, (ny+1)/2 - 1
1272	f_vec_y(mm,2*j+1) = ar(ny+1-j,i,k)
1273	ENDDO
1274	ENDDO
1275	ENDDO
1276
1277	ENDIF
1278
1279	f_vec_y(:,1) = 0.0_wp
1280	f_vec_y(:,ny+2) = 0.0_wp
1281
1282	ALLOCATE( work_vec((nxr_y_l-nxl_y_l+1)*(nzt_y-nzb_y+1),ny+2) )
1283	CALL fft991cy_vec( f_vec_y, work_vec, trigs_y, ifax_y, ny+1, 1 )
1284	DEALLOCATE( work_vec )
1285
1286	mm = 1
1287	DO k = nzb_y, nzt_y
1288	DO i = nxl_y_l, nxr_y_l
1289	ar(0:ny,i,k) = f_vec_y(mm,0:ny)
1290	mm = mm+1
1291	ENDDO
1292	ENDDO
1293
1294	DEALLOCATE( f_vec_y )
1295
1296	ENDIF
1297
1298	ENDIF
1299
1300	ELSEIF ( fft_method == 'fftw' ) THEN
1301
1302	#if defined( __fftw )
1303	IF ( forward_fft ) THEN
1304
1305	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1306	!$OMP DO
1307	DO k = nzb_y, nzt_y
1308	DO i = nxl_y_l, nxr_y_l
1309
1310	y_in(0:ny) = ar(0:ny,i,k)
1311	CALL FFTW_EXECUTE_DFT_R2C( plan_yf, y_in, y_out )
1312
1313	DO j = 0, (ny+1)/2
1314	ar_tr(j,i,k) = REAL( y_out(j), KIND=wp ) / (ny+1)
1315	ENDDO
1316	DO j = 1, (ny+1)/2 - 1
1317	ar_tr(ny+1-j,i,k) = AIMAG( y_out(j) ) / (ny+1)
1318	ENDDO
1319
1320	ENDDO
1321	ENDDO
1322	!$OMP END PARALLEL
1323
1324	ELSE
1325
1326	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1327	!$OMP DO
1328	DO k = nzb_y, nzt_y
1329	DO i = nxl_y_l, nxr_y_l
1330
1331	y_out(0) = CMPLX( ar_tr(0,i,k), 0.0_wp, KIND=wp )
1332	DO j = 1, (ny+1)/2 - 1
1333	y_out(j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k), KIND=wp )
1334	ENDDO
1335	y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp, KIND=wp )
1336
1337	CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
1338	ar(0:ny,i,k) = y_in(0:ny)
1339
1340	ENDDO
1341	ENDDO
1342	!$OMP END PARALLEL
1343
1344	ENDIF
1345	#endif
1346
1347	ELSEIF ( fft_method == 'system-specific' ) THEN
1348
1349	#if defined( __ibm )
1350	IF ( forward_fft) THEN
1351
1352	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1353	!$OMP DO
1354	DO k = nzb_y, nzt_y
1355	DO i = nxl_y_l, nxr_y_l
1356
1357	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, auy2, nau2 )
1358
1359	DO j = 0, (ny+1)/2
1360	ar_tr(j,i,k) = work(2*j)
1361	ENDDO
1362	DO j = 1, (ny+1)/2 - 1
1363	ar_tr(ny+1-j,i,k) = work(2*j+1)
1364	ENDDO
1365
1366	ENDDO
1367	ENDDO
1368	!$OMP END PARALLEL
1369
1370	ELSE
1371
1372	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1373	!$OMP DO
1374	DO k = nzb_y, nzt_y
1375	DO i = nxl_y_l, nxr_y_l
1376
1377	DO j = 0, (ny+1)/2
1378	work(2*j) = ar_tr(j,i,k)
1379	ENDDO
1380	DO j = 1, (ny+1)/2 - 1
1381	work(2*j+1) = ar_tr(ny+1-j,i,k)
1382	ENDDO
1383	work(1) = 0.0_wp
1384	work(ny+2) = 0.0_wp
1385
1386	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, auy4, nau2 )
1387
1388	DO j = 0, ny
1389	ar(j,i,k) = work(j)
1390	ENDDO
1391
1392	ENDDO
1393	ENDDO
1394	!$OMP END PARALLEL
1395
1396	ENDIF
1397	#elif defined( __nec_fft )
1398	IF ( forward_fft ) THEN
1399
1400	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1401	!$OMP DO
1402	DO k = nzb_y, nzt_y
1403	DO i = nxl_y_l, nxr_y_l
1404
1405	work(0:ny) = ar(0:ny,i,k)
1406
1407	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
1408
1409	DO j = 0, (ny+1)/2
1410	ar_tr(j,i,k) = work(2*j)
1411	ENDDO
1412	DO j = 1, (ny+1)/2 - 1
1413	ar_tr(ny+1-j,i,k) = work(2*j+1)
1414	ENDDO
1415
1416	ENDDO
1417	ENDDO
1418	!$END OMP PARALLEL
1419
1420	ELSE
1421
1422	!$OMP PARALLEL PRIVATE ( work, i, j, k )
1423	!$OMP DO
1424	DO k = nzb_y, nzt_y
1425	DO i = nxl_y_l, nxr_y_l
1426
1427	DO j = 0, (ny+1)/2
1428	work(2*j) = ar_tr(j,i,k)
1429	ENDDO
1430	DO j = 1, (ny+1)/2 - 1
1431	work(2*j+1) = ar_tr(ny+1-j,i,k)
1432	ENDDO
1433	work(1) = 0.0_wp
1434	work(ny+2) = 0.0_wp
1435
1436	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1437
1438	ar(0:ny,i,k) = work(0:ny)
1439
1440	ENDDO
1441	ENDDO
1442	!$OMP END PARALLEL
1443
1444	ENDIF
1445	#elif defined( __cuda_fft )
1446
1447	IF ( forward_fft ) THEN
1448
1449	!$ACC HOST_DATA USE_DEVICE(ar, ar_tmp)
1450	CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
1451	!$ACC END HOST_DATA
1452
1453	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i,j,k) &
1454	!$ACC PRESENT(ar, ar_tmp)
1455	DO k = nzb_y, nzt_y
1456	DO i = nxl_y, nxr_y
1457
1458	DO j = 0, (ny+1)/2
1459	ar(j,i,k) = REAL( ar_tmp(j,i,k), KIND=wp ) * dny
1460	ENDDO
1461
1462	DO j = 1, (ny+1)/2 - 1
1463	ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
1464	ENDDO
1465
1466	ENDDO
1467	ENDDO
1468
1469	ELSE
1470
1471	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i,j,k) &
1472	!$ACC PRESENT(ar, ar_tmp)
1473	DO k = nzb_y, nzt_y
1474	DO i = nxl_y, nxr_y
1475
1476	ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0_wp, KIND=wp )
1477
1478	DO j = 1, (ny+1)/2 - 1
1479	ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k), KIND=wp )
1480	ENDDO
1481	ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0_wp, KIND=wp )
1482
1483	ENDDO
1484	ENDDO
1485
1486	!$ACC HOST_DATA USE_DEVICE(ar, ar_tmp)
1487	CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
1488	!$ACC END HOST_DATA
1489
1490	ENDIF
1491
1492	#endif
1493
1494	ENDIF
1495
1496	END SUBROUTINE fft_y
1497
1498	!--------------------------------------------------------------------------------------------------!
1499	! Description:
1500	! ------------
1501	!> Fourier-transformation along y-direction.
1502	!> Version for 1D-decomposition.
1503	!> It uses internal algorithms (Singleton or Temperton) or system-specific routines, if they are
1504	!> available.
1505	!--------------------------------------------------------------------------------------------------!
1506
1507	SUBROUTINE fft_y_1d( ar, direction )
1508
1509
1510	IMPLICIT NONE
1511
1512	CHARACTER (LEN=*) :: direction
1513
1514	INTEGER(iwp) :: j !<
1515	INTEGER(iwp) :: jshape(1) !<
1516
1517	LOGICAL :: forward_fft !<
1518
1519	REAL(wp), DIMENSION(0:ny) :: ar !<
1520	REAL(wp), DIMENSION(0:ny+2) :: work !<
1521	REAL(wp), DIMENSION(ny+2) :: work1 !<
1522
1523	COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: cwork !<
1524
1525	#if defined( __ibm )
1526	REAL(wp), DIMENSION(nau2) :: auy2 !<
1527	REAL(wp), DIMENSION(nau2) :: auy4 !<
1528	#elif defined( __nec_fft )
1529	REAL(wp), DIMENSION(6*(ny+1)) :: work2 !<
1530	#endif
1531
1532	IF ( direction == 'forward' ) THEN
1533	forward_fft = .TRUE.
1534	ELSE
1535	forward_fft = .FALSE.
1536	ENDIF
1537
1538	IF ( fft_method == 'singleton-algorithm' ) THEN
1539
1540	!
1541	!-- Performing the fft with singleton's software works on every system, since it is part of
1542	!-- the model.
1543	ALLOCATE( cwork(0:ny) )
1544
1545	IF ( forward_fft ) THEN
1546
1547	DO j = 0, ny
1548	cwork(j) = CMPLX( ar(j), KIND=wp )
1549	ENDDO
1550
1551	jshape = SHAPE( cwork )
1552	CALL FFTN( cwork, jshape )
1553
1554	DO j = 0, (ny+1)/2
1555	ar(j) = REAL( cwork(j), KIND=wp )
1556	ENDDO
1557	DO j = 1, (ny+1)/2 - 1
1558	ar(ny+1-j) = -AIMAG( cwork(j) )
1559	ENDDO
1560
1561	ELSE
1562
1563	cwork(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
1564	DO j = 1, (ny+1)/2 - 1
1565	cwork(j) = CMPLX( ar(j), -ar(ny+1-j), KIND=wp )
1566	cwork(ny+1-j) = CMPLX( ar(j), ar(ny+1-j), KIND=wp )
1567	ENDDO
1568	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp, KIND=wp )
1569
1570	jshape = SHAPE( cwork )
1571	CALL FFTN( cwork, jshape, inv = .TRUE. )
1572
1573	DO j = 0, ny
1574	ar(j) = REAL( cwork(j), KIND=wp )
1575	ENDDO
1576
1577	ENDIF
1578
1579	DEALLOCATE( cwork )
1580
1581	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
1582
1583	!
1584	!-- Performing the fft with Temperton's software works on every system, since it is part of
1585	!-- the model.
1586	IF ( forward_fft ) THEN
1587
1588	work(0:ny) = ar
1589	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
1590
1591	DO j = 0, (ny+1)/2
1592	ar(j) = work(2*j)
1593	ENDDO
1594	DO j = 1, (ny+1)/2 - 1
1595	ar(ny+1-j) = work(2*j+1)
1596	ENDDO
1597
1598	ELSE
1599
1600	DO j = 0, (ny+1)/2
1601	work(2*j) = ar(j)
1602	ENDDO
1603	DO j = 1, (ny+1)/2 - 1
1604	work(2*j+1) = ar(ny+1-j)
1605	ENDDO
1606	work(1) = 0.0_wp
1607	work(ny+2) = 0.0_wp
1608
1609	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
1610	ar = work(0:ny)
1611
1612	ENDIF
1613
1614	ELSEIF ( fft_method == 'fftw' ) THEN
1615
1616	#if defined( __fftw )
1617	IF ( forward_fft ) THEN
1618
1619	y_in(0:ny) = ar(0:ny)
1620	CALL FFTW_EXECUTE_DFT_R2C( plan_yf, y_in, y_out )
1621
1622	DO j = 0, (ny+1)/2
1623	ar(j) = REAL( y_out(j), KIND=wp ) / (ny+1)
1624	ENDDO
1625	DO j = 1, (ny+1)/2 - 1
1626	ar(ny+1-j) = AIMAG( y_out(j) ) / (ny+1)
1627	ENDDO
1628
1629	ELSE
1630
1631	y_out(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
1632	DO j = 1, (ny+1)/2 - 1
1633	y_out(j) = CMPLX( ar(j), ar(ny+1-j), KIND=wp )
1634	ENDDO
1635	y_out((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp, KIND=wp )
1636
1637	CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
1638	ar(0:ny) = y_in(0:ny)
1639
1640	ENDIF
1641	#endif
1642
1643	ELSEIF ( fft_method == 'system-specific' ) THEN
1644
1645	#if defined( __ibm )
1646	IF ( forward_fft ) THEN
1647
1648	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, auy2, nau2 )
1649
1650	DO j = 0, (ny+1)/2
1651	ar(j) = work(2*j)
1652	ENDDO
1653	DO j = 1, (ny+1)/2 - 1
1654	ar(ny+1-j) = work(2*j+1)
1655	ENDDO
1656
1657	ELSE
1658
1659	DO j = 0, (ny+1)/2
1660	work(2*j) = ar(j)
1661	ENDDO
1662	DO j = 1, (ny+1)/2 - 1
1663	work(2*j+1) = ar(ny+1-j)
1664	ENDDO
1665	work(1) = 0.0_wp
1666	work(ny+2) = 0.0_wp
1667
1668	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, auy4, nau2 )
1669
1670	DO j = 0, ny
1671	ar(j) = work(j)
1672	ENDDO
1673
1674	ENDIF
1675	#elif defined( __nec_fft )
1676	IF ( forward_fft ) THEN
1677
1678	work(0:ny) = ar(0:ny)
1679
1680	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
1681
1682	DO j = 0, (ny+1)/2
1683	ar(j) = work(2*j)
1684	ENDDO
1685	DO j = 1, (ny+1)/2 - 1
1686	ar(ny+1-j) = work(2*j+1)
1687	ENDDO
1688
1689	ELSE
1690
1691	DO j = 0, (ny+1)/2
1692	work(2*j) = ar(j)
1693	ENDDO
1694	DO j = 1, (ny+1)/2 - 1
1695	work(2*j+1) = ar(ny+1-j)
1696	ENDDO
1697	work(1) = 0.0_wp
1698	work(ny+2) = 0.0_wp
1699
1700	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1701
1702	ar(0:ny) = work(0:ny)
1703
1704	ENDIF
1705	#endif
1706
1707	ENDIF
1708
1709	END SUBROUTINE fft_y_1d
1710
1711	!--------------------------------------------------------------------------------------------------!
1712	! Description:
1713	! ------------
1714	!> Fourier-transformation along x-direction.
1715	!> Version for 1d domain decomposition,
1716	!> using multiple 1D FFT from Math Keisan on NEC or Temperton-algorithm
1717	!> (no singleton-algorithm on NEC because it does not vectorize).
1718	!--------------------------------------------------------------------------------------------------!
1719
1720	SUBROUTINE fft_x_m( ar, direction )
1721
1722
1723	IMPLICIT NONE
1724
1725	CHARACTER (LEN=*) :: direction !<
1726
1727	INTEGER(iwp) :: i !<
1728	INTEGER(iwp) :: k !<
1729	INTEGER(iwp) :: siza !<
1730	#if defined( __nec_fft )
1731	INTEGER(iwp) :: sizw
1732	#endif
1733
1734	REAL(wp), DIMENSION(0:nx,nz) :: ar !<
1735	REAL(wp), DIMENSION(0:nx+3,nz+1) :: ai !<
1736	REAL(wp), DIMENSION(6*(nx+4),nz+1) :: work1 !<
1737
1738	#if defined( __nec_fft )
1739	COMPLEX(wp), DIMENSION(:,:), ALLOCATABLE :: work
1740	#endif
1741
1742	IF ( fft_method == 'temperton-algorithm' ) THEN
1743
1744	siza = SIZE( ai, 1 )
1745
1746	IF ( direction == 'forward') THEN
1747
1748	ai(0:nx,1:nz) = ar(0:nx,1:nz)
1749	ai(nx+1:,:) = 0.0_wp
1750
1751	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
1752
1753	DO k = 1, nz
1754	DO i = 0, (nx+1)/2
1755	ar(i,k) = ai(2*i,k)
1756	ENDDO
1757	DO i = 1, (nx+1)/2 - 1
1758	ar(nx+1-i,k) = ai(2*i+1,k)
1759	ENDDO
1760	ENDDO
1761
1762	ELSE
1763
1764	DO k = 1, nz
1765	DO i = 0, (nx+1)/2
1766	ai(2*i,k) = ar(i,k)
1767	ENDDO
1768	DO i = 1, (nx+1)/2 - 1
1769	ai(2*i+1,k) = ar(nx+1-i,k)
1770	ENDDO
1771	ai(1,k) = 0.0_wp
1772	ai(nx+2,k) = 0.0_wp
1773	ENDDO
1774
1775	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
1776
1777	ar(0:nx,1:nz) = ai(0:nx,1:nz)
1778
1779	ENDIF
1780
1781	ELSEIF ( fft_method == 'system-specific' ) THEN
1782
1783	#if defined( __nec_fft )
1784	ALLOCATE( work((nx+4)/2+1,nz+1) )
1785	siza = SIZE( ai, 1 )
1786	sizw = SIZE( work, 1 )
1787
1788	IF ( direction == 'forward') THEN
1789
1790	!
1791	!-- Tables are initialized once more. This call should not be necessary, but otherwise
1792	!-- program aborts in asymmetric case.
1793	CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, trig_xf, work1, 0 )
1794
1795	ai(0:nx,1:nz) = ar(0:nx,1:nz)
1796	IF ( nz1 > nz ) THEN
1797	ai(:,nz1) = 0.0_wp
1798	ENDIF
1799
1800	CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw, trig_xf, work1, 0 )
1801
1802	DO k = 1, nz
1803	DO i = 0, (nx+1)/2
1804	ar(i,k) = REAL( work(i+1,k), KIND=wp )
1805	ENDDO
1806	DO i = 1, (nx+1)/2 - 1
1807	ar(nx+1-i,k) = AIMAG( work(i+1,k) )
1808	ENDDO
1809	ENDDO
1810
1811	ELSE
1812
1813	!
1814	!-- Tables are initialized once more. This call should not be
1815	!-- necessary, but otherwise program aborts in asymmetric case
1816	CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, trig_xb, work1, 0 )
1817
1818	IF ( nz1 > nz ) THEN
1819	work(:,nz1) = 0.0_wp
1820	ENDIF
1821	DO k = 1, nz
1822	work(1,k) = CMPLX( ar(0,k), 0.0_wp, KIND=wp )
1823	DO i = 1, (nx+1)/2 - 1
1824	work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k), KIND=wp )
1825	ENDDO
1826	work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0_wp, KIND=wp )
1827	ENDDO
1828
1829	CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, trig_xb, work1, 0 )
1830
1831	ar(0:nx,1:nz) = ai(0:nx,1:nz)
1832
1833	ENDIF
1834
1835	DEALLOCATE( work )
1836	#endif
1837
1838	ENDIF
1839
1840	END SUBROUTINE fft_x_m
1841
1842	!--------------------------------------------------------------------------------------------------!
1843	! Description:
1844	! ------------
1845	!> Fourier-transformation along y-direction.
1846	!> Version for 1d domain decomposition,
1847	!> using multiple 1D FFT from Math Keisan on NEC or Temperton-algorithm
1848	!> (no singleton-algorithm on NEC because it does not vectorize).
1849	!--------------------------------------------------------------------------------------------------!
1850
1851	SUBROUTINE fft_y_m( ar, ny1, direction )
1852
1853
1854	IMPLICIT NONE
1855
1856	CHARACTER (LEN=*) :: direction !<
1857
1858	INTEGER(iwp) :: j !<
1859	INTEGER(iwp) :: k !<
1860	INTEGER(iwp) :: ny1 !<
1861	INTEGER(iwp) :: siza !<
1862	#if defined( __nec_fft )
1863	INTEGER(iwp) :: sizw
1864	#endif
1865
1866	REAL(wp), DIMENSION(0:ny1,nz) :: ar !<
1867	REAL(wp), DIMENSION(0:ny+3,nz+1) :: ai !<
1868	REAL(wp), DIMENSION(6*(ny+4),nz+1) :: work1 !<
1869
1870	#if defined( __nec_fft )
1871	COMPLEX(wp), DIMENSION(:,:), ALLOCATABLE :: work
1872	#endif
1873
1874
1875	IF ( fft_method == 'temperton-algorithm' ) THEN
1876
1877	siza = SIZE( ai, 1 )
1878
1879	IF ( direction == 'forward') THEN
1880
1881	ai(0:ny,1:nz) = ar(0:ny,1:nz)
1882	ai(ny+1:,:) = 0.0_wp
1883
1884	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
1885
1886	DO k = 1, nz
1887	DO j = 0, (ny+1)/2
1888	ar(j,k) = ai(2*j,k)
1889	ENDDO
1890	DO j = 1, (ny+1)/2 - 1
1891	ar(ny+1-j,k) = ai(2*j+1,k)
1892	ENDDO
1893	ENDDO
1894
1895	ELSE
1896
1897	DO k = 1, nz
1898	DO j = 0, (ny+1)/2
1899	ai(2*j,k) = ar(j,k)
1900	ENDDO
1901	DO j = 1, (ny+1)/2 - 1
1902	ai(2*j+1,k) = ar(ny+1-j,k)
1903	ENDDO
1904	ai(1,k) = 0.0_wp
1905	ai(ny+2,k) = 0.0_wp
1906	ENDDO
1907
1908	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
1909
1910	ar(0:ny,1:nz) = ai(0:ny,1:nz)
1911
1912	ENDIF
1913
1914	ELSEIF ( fft_method == 'system-specific' ) THEN
1915
1916	#if defined( __nec_fft )
1917	ALLOCATE( work((ny+4)/2+1,nz+1) )
1918	siza = SIZE( ai, 1 )
1919	sizw = SIZE( work, 1 )
1920
1921	IF ( direction == 'forward') THEN
1922
1923	!
1924	!-- Tables are initialized once more. This call should not be necessary, but otherwise
1925	!-- program aborts in asymmetric case.
1926	CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, trig_yf, work1, 0 )
1927
1928	ai(0:ny,1:nz) = ar(0:ny,1:nz)
1929	IF ( nz1 > nz ) THEN
1930	ai(:,nz1) = 0.0_wp
1931	ENDIF
1932
1933	CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, trig_yf, work1, 0 )
1934
1935	DO k = 1, nz
1936	DO j = 0, (ny+1)/2
1937	ar(j,k) = REAL( work(j+1,k), KIND=wp )
1938	ENDDO
1939	DO j = 1, (ny+1)/2 - 1
1940	ar(ny+1-j,k) = AIMAG( work(j+1,k) )
1941	ENDDO
1942	ENDDO
1943
1944	ELSE
1945
1946	!
1947	!-- Tables are initialized once more. This call should not be necessary, but otherwise
1948	!-- program aborts in asymmetric case.
1949	CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, trig_yb, work1, 0 )
1950
1951	IF ( nz1 > nz ) THEN
1952	work(:,nz1) = 0.0_wp
1953	ENDIF
1954	DO k = 1, nz
1955	work(1,k) = CMPLX( ar(0,k), 0.0_wp, KIND=wp )
1956	DO j = 1, (ny+1)/2 - 1
1957	work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k), KIND=wp )
1958	ENDDO
1959	work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0_wp, KIND=wp )
1960	ENDDO
1961
1962	CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, trig_yb, work1, 0 )
1963
1964	ar(0:ny,1:nz) = ai(0:ny,1:nz)
1965
1966	ENDIF
1967
1968	DEALLOCATE( work )
1969	#endif
1970
1971	ENDIF
1972
1973	END SUBROUTINE fft_y_m
1974
1975
1976	END MODULE fft_xy

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |