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

Last change on this file since 1354 was 1354, checked in by heinze, 8 years ago

last commit documented

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