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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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