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

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

last commit documented

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