Home

Context Navigation

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

Last change on this file since 3039 was 2718, checked in by maronga, 7 years ago
deleting of deprecated files; headers updated where needed
Property svn:keywords set to `Id`
File size: 53.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |