Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |