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

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

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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