Home

Context Navigation

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

Last change on this file since 3774 was 3655, checked in by knoop, 6 years ago
Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update
Property svn:keywords set to `Id`
File size: 55.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |