Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |