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

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

last commit documented

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