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

Last change on this file since 1482 was 1482, checked in by raasch, 7 years ago

adjustments for using CUDA-aware MPI

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