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

Last change on this file since 1434 was 1403, checked in by raasch, 10 years ago

last commit documented

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