Home

Context Navigation

fft_xy.f90 @ 1106

Last change on this file since 1106 was 1106, checked in by raasch, 11 years ago

New:
---

Porting of FFT-solver for serial runs to GPU using CUDA FFT,
preprocessor lines in transpose routines rearranged, so that routines can also
be used in serial (non-parallel) mode,
transpositions also carried out in serial mode, routines fftx, fftxp replaced
by calls of fft_x, fft_x replaced by fft_x_1d in the 1D-decomposition routines
(Makefile, Makefile_check, fft_xy, poisfft, poisfft_hybrid, transpose, new: cuda_fft_interfaces)

--stdin argument for mpiexec on lckyuh, -y and -Y settings output to header (mrun)

Changed:

Module array_kind renamed precision_kind
(check_open, data_output_3d, fft_xy, modules, user_data_output_3d)

some format changes for coupled atmosphere-ocean runs (header)
small changes in code formatting (microphysics, prognostic_equations)

Errors:

bugfix: default value (0) assigned to coupling_start_time (modules)
bugfix: initial time for preruns of coupled runs is output as -coupling_start_time (data_output_profiles)

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

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

Context Navigation

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

Download in other formats: