Home

Context Navigation

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

Last change on this file since 2000 was 2000, checked in by knoop, 8 years ago
Forced header and separation lines into 80 columns
Property svn:keywords set to `Id`
File size: 57.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |