Home

Context Navigation

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

Last change on this file since 2416 was 2300, checked in by raasch, 7 years ago
NEC related code partly removed, host variable partly removed, host specific code completely removed, default values for host, loop_optimization and termination time_needed changed
Property svn:keywords set to `Id`
File size: 53.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |