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

Last change on this file since 1818 was 1818, checked in by maronga, 5 years ago

last commit documented / copyright update

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