Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |