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

Last change on this file since 1375 was 1375, checked in by raasch, 7 years ago

last commit documented

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