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

Last change on this file since 1267 was 1258, checked in by raasch, 11 years ago

last commit documented

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