Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |