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

Last change on this file since 1310 was 1310, checked in by raasch, 10 years ago

update of GPL copyright

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