Home

Context Navigation

source: palm/trunk/SOURCE/fft_xy.f90 @ 1196

Last change on this file since 1196 was 1167, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 46.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |