Home

Context Navigation

fft_xy_mod.f90 @ 4869

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

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

Context Navigation

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

Download in other formats: