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

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

openMP-bugfix for fftw: some arrays defined as threadprivate

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