Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |