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

Last change on this file since 1762 was 1750, checked in by raasch, 8 years ago

last commit documented

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