Home

Context Navigation

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

Last change on this file since 2131 was 2119, checked in by raasch, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 52.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |