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

Last change on this file since 1713 was 1683, checked in by knoop, 9 years ago

last commit documented

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