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

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

last commit documented

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