Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |