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

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

last commit documented

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