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

Last change on this file since 1850 was 1850, checked in by maronga, 5 years ago

added _mod string to several filenames to meet the naming convection for modules

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