source: palm/tags/release-4.0/SOURCE/fft_xy.f90 @ 4239

Last change on this file since 4239 was 1483, checked in by raasch, 9 years ago

last commit documented

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