Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |