source: palm/trunk/SOURCE/fft_xy_mod.f90 @ 2036

Last change on this file since 2036 was 2001, checked in by knoop, 8 years ago

last commit documented

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