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

Last change on this file since 1323 was 1323, checked in by raasch, 8 years ago

last commit documented

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