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
RevLine 
[1]1 MODULE fft_xy
2
[1036]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!
[254]20! Current revisions:
[1]21! -----------------
22!
[1258]23!
[1]24! Former revisions:
25! -----------------
[3]26! $Id: fft_xy.f90 1258 2013-11-08 16:09:09Z heinze $
[392]27!
[1258]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!
[1220]32! 1219 2013-08-30 09:33:18Z heinze
33! bugfix: use own branch for fftw
34!
[1217]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!
[1211]40! 1210 2013-08-14 10:58:20Z raasch
41! fftw added
42!
[1167]43! 1166 2013-05-24 13:55:44Z raasch
44! C_DOUBLE/COMPLEX reset to dpk
45!
[1154]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!
[1112]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!
[1107]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!
[1093]59! 1092 2013-02-02 11:24:22Z raasch
60! variable sizw declared for NEC case only
61!
[1037]62! 1036 2012-10-22 13:43:42Z raasch
63! code put under GPL (PALM 3.9)
64!
[392]65! 274 2009-03-26 15:11:21Z heinze
66! Output of messages replaced by message handling routine.
67!
68! Feb. 2007
[3]69! RCS Log replace by Id keyword, revision history cleaned up
70!
[1]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
[1153]97#if defined( __cuda_fft )
98    USE ISO_C_BINDING
[1210]99#elif defined( __fftw )
100    USE, INTRINSIC ::  ISO_C_BINDING
[1153]101#endif
[1106]102    USE precision_kind
[1]103    USE singleton
104    USE temperton_fft
[1106]105    USE transpose_indices
[1]106
107    IMPLICIT NONE
108
109    PRIVATE
[1106]110    PUBLIC fft_x, fft_x_1d, fft_y, fft_y_1d, fft_init, fft_x_m, fft_y_m
[1]111
112    INTEGER, DIMENSION(:), ALLOCATABLE, SAVE ::  ifax_x, ifax_y
113
114    LOGICAL, SAVE                            ::  init_fft = .FALSE.
115
[1106]116    REAL, SAVE ::  dnx, dny, sqr_dnx, sqr_dny
[1]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
[1106]129#elif defined( __cuda_fft )
[1153]130    INTEGER(C_INT), SAVE ::  plan_xf, plan_xi, plan_yf, plan_yi
131    INTEGER, SAVE        ::  total_points_x_transpo, total_points_y_transpo
[1219]132#endif
133
134#if defined( __fftw )
[1210]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
[1]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
[1106]152    INTERFACE fft_x_1d
153       MODULE PROCEDURE fft_x_1d
154    END INTERFACE fft_x_1d
155
[1]156    INTERFACE fft_y
157       MODULE PROCEDURE fft_y
158    END INTERFACE fft_y
159
[1106]160    INTERFACE fft_y_1d
161       MODULE PROCEDURE fft_y_1d
162    END INTERFACE fft_y_1d
163
[1]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
[1106]177       USE cuda_fft_interfaces
178
[1]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
[1106]205          dnx = 1.0 / ( nx + 1.0 )
206          dny = 1.0 / ( ny + 1.0 )
207          sqr_dnx = SQRT( dnx )
208          sqr_dny = SQRT( dny )
[1]209#if defined( __ibm )  &&  ! defined( __ibmy_special )
210!
211!--       Initialize tables for fft along x
[1106]212          CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1,  1, sqr_dnx, aux1, nau1, &
[1]213                      aux2, nau2 )
[1106]214          CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
[1]215                      aux4, nau2 )
216!
217!--       Initialize tables for fft along y
[1106]218          CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1,  1, sqr_dny, auy1, nau1, &
[1]219                      auy2, nau2 )
[1106]220          CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
[1]221                      auy4, nau2 )
222#elif defined( __nec )
[254]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 )
[1]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))
[1106]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, &
[1]240                       trig_xf, workx, 0 )
[1106]241          CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
[1]242                       trig_xb, workx, 0 )
243!
244!--       Initialize tables for fft along y (non-vector and vector case (M))
[1106]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, &
[1]248                       trig_yf, worky, 0 )
[1106]249          CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
[1]250                       trig_yb, worky, 0 )
[1106]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)
[1111]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) )
[1]258#else
[254]259          message_string = 'no system-specific fft-call available'
260          CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
[1]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
[1210]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
[1]288       ELSEIF ( fft_method == 'singleton-algorithm' )  THEN
289
290          CONTINUE
291
292       ELSE
293
[254]294          message_string = 'fft method "' // TRIM( fft_method) // &
295                           '" not available'
296          CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
[1]297       ENDIF
298
299    END SUBROUTINE fft_init
300
301
[1216]302    SUBROUTINE fft_x( ar, direction, ar_2d )
[1]303
304!----------------------------------------------------------------------!
305!                                 fft_x                                !
306!                                                                      !
307!               Fourier-transformation along x-direction               !
[1106]308!                     Version for 2D-decomposition                     !
[1]309!                                                                      !
310!      fft_x uses internal algorithms (Singleton or Temperton) or      !
311!           system-specific routines, if they are available            !
312!----------------------------------------------------------------------!
313
[1106]314       USE cuda_fft_interfaces
[1153]315#if defined( __cuda_fft )
316       USE ISO_C_BINDING
317#endif
[1106]318
[1]319       IMPLICIT NONE
320
321       CHARACTER (LEN=*) ::  direction
[1111]322       INTEGER ::  i, ishape(1), j, k
[1106]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 )
[1257]334       COMPLEX(dpk), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) ::  ar_tmp
[1111]335       !$acc declare create( ar_tmp )
[1106]336#endif
[1216]337       REAL, DIMENSION(0:nx,nys_x:nyn_x), OPTIONAL   ::  ar_2d
[1106]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
[1210]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
[1216]471                   IF ( PRESENT( ar_2d ) )  THEN
[1210]472
[1216]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
[1210]491                ENDDO
492             ENDDO
493             !$OMP END PARALLEL
494
[1216]495          ELSE
[1210]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
[1216]501                   IF ( PRESENT( ar_2d ) )  THEN
[1210]502
[1216]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
[1210]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
[1216]526          ENDIF
[1210]527#endif
528
[1106]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
[1111]636             !$acc data present( ar )
637             CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
[1106]638
[1111]639             !$acc kernels
[1106]640             DO  k = nzb_x, nzt_x
641                DO  j = nys_x, nyn_x
642
643                   DO  i = 0, (nx+1)/2
[1111]644                      ar(i,j,k)      = REAL( ar_tmp(i,j,k) )  * dnx
[1106]645                   ENDDO
646
647                   DO  i = 1, (nx+1)/2 - 1
[1111]648                      ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
[1106]649                   ENDDO
650
651                ENDDO
652             ENDDO
[1111]653             !$acc end kernels
654             !$acc end data
[1106]655
656          ELSE
657
[1111]658             !$acc data present( ar )
659             !$acc kernels
[1106]660             DO  k = nzb_x, nzt_x
661                DO  j = nys_x, nyn_x
662
[1111]663                   ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
[1106]664
665                   DO  i = 1, (nx+1)/2 - 1
[1111]666                      ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
[1106]667                   ENDDO
[1111]668                   ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
[1106]669
670                ENDDO
671             ENDDO
[1111]672             !$acc end kernels
[1106]673
[1111]674             CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
675             !$acc end data
[1106]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
[1]709       INTEGER ::  i, ishape(1)
710
[1106]711       LOGICAL ::  forward_fft
712
[1]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
[1106]723       IF ( direction == 'forward' )  THEN
724          forward_fft = .TRUE.
725       ELSE
726          forward_fft = .FALSE.
727       ENDIF
728
[1]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     
[1106]736          IF ( forward_fft )   then
[1]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
[1106]775          IF ( forward_fft )  THEN
[1]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
[1216]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
[1]832       ELSEIF ( fft_method == 'system-specific' )  THEN
833
834#if defined( __ibm )  &&  ! defined( __ibmy_special )
[1106]835          IF ( forward_fft )  THEN
[1]836
[1106]837             CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
[1]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
[1106]858             CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
[1]859                         aux4, nau2 )
860
861             DO  i = 0, nx
862                ar(i) = work(i)
863             ENDDO
864
865          ENDIF
866#elif defined( __nec )
[1106]867          IF ( forward_fft )  THEN
[1]868
869             work(0:nx) = ar(0:nx)
870
[1106]871             CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
872     
[1]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
[1106]891             CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
[1]892
893             ar(0:nx) = work(0:nx)
894
895          ENDIF
896#else
[254]897          message_string = 'no system-specific fft-call available'
[1106]898          CALL message( 'fft_x_1d', 'PA0188', 1, 2, 0, 6, 0 )
[1]899#endif
900       ELSE
[274]901          message_string = 'fft method "' // TRIM( fft_method) // &
902                           '" not available'
[1106]903          CALL message( 'fft_x_1d', 'PA0189', 1, 2, 0, 6, 0 )
[1]904
905       ENDIF
906
[1106]907    END SUBROUTINE fft_x_1d
[1]908
[1216]909    SUBROUTINE fft_y( ar, direction, ar_tr, nxl_y_bound, nxr_y_bound, nxl_y_l, &
910                      nxr_y_l )
[1]911
912!----------------------------------------------------------------------!
913!                                 fft_y                                !
914!                                                                      !
915!               Fourier-transformation along y-direction               !
[1106]916!                     Version for 2D-decomposition                     !
[1]917!                                                                      !
918!      fft_y uses internal algorithms (Singleton or Temperton) or      !
919!           system-specific routines, if they are available            !
[1216]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.                 !
[1]935!----------------------------------------------------------------------!
936
[1106]937       USE cuda_fft_interfaces
[1153]938#if defined( __cuda_fft )
939       USE ISO_C_BINDING
940#endif
[1106]941
[1]942       IMPLICIT NONE
943
944       CHARACTER (LEN=*) ::  direction
[1111]945       INTEGER ::  i, j, jshape(1), k
[1216]946       INTEGER ::  nxl_y_bound, nxl_y_l, nxr_y_bound, nxr_y_l
[1106]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 )
[1257]958       COMPLEX(dpk), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) ::  ar_tmp
[1111]959       !$acc declare create( ar_tmp )
[1106]960#endif
[1216]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
[1106]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
[1216]982                DO  i = nxl_y_l, nxr_y_l
[1106]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
[1216]992                      ar_tr(j,i,k) = REAL( cwork(j) )
[1106]993                   ENDDO
994                   DO  j = 1, (ny+1)/2 - 1
[1216]995                      ar_tr(ny+1-j,i,k) = -AIMAG( cwork(j) )
[1106]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
[1216]1007                DO  i = nxl_y_l, nxr_y_l
[1106]1008
[1216]1009                   cwork(0) = CMPLX( ar_tr(0,i,k), 0.0 )
[1106]1010                   DO  j = 1, (ny+1)/2 - 1
[1216]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) )
[1106]1013                   ENDDO
[1216]1014                   cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 )
[1106]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
[1216]1041                DO  i = nxl_y_l, nxr_y_l
[1106]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
[1216]1047                      ar_tr(j,i,k) = work(2*j)
[1106]1048                   ENDDO
1049                   DO  j = 1, (ny+1)/2 - 1
[1216]1050                      ar_tr(ny+1-j,i,k) = work(2*j+1)
[1106]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
[1216]1062                DO  i = nxl_y_l, nxr_y_l
[1106]1063
1064                   DO  j = 0, (ny+1)/2
[1216]1065                      work(2*j) = ar_tr(j,i,k)
[1106]1066                   ENDDO
1067                   DO  j = 1, (ny+1)/2 - 1
[1216]1068                      work(2*j+1) = ar_tr(ny+1-j,i,k)
[1106]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
[1210]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
[1216]1090                DO  i = nxl_y_l, nxr_y_l
[1210]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
[1216]1096                      ar_tr(j,i,k) = REAL( y_out(j) ) / (ny+1)
[1210]1097                   ENDDO
1098                   DO  j = 1, (ny+1)/2 - 1
[1216]1099                      ar_tr(ny+1-j,i,k) = AIMAG( y_out(j) ) / (ny+1)
[1210]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
[1216]1111                DO  i = nxl_y_l, nxr_y_l
[1210]1112
[1216]1113                   y_out(0) = CMPLX( ar_tr(0,i,k), 0.0 )
[1210]1114                   DO  j = 1, (ny+1)/2 - 1
[1216]1115                      y_out(j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k) )
[1210]1116                   ENDDO
[1216]1117                   y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 )
[1210]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
[1106]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
[1216]1137                DO  i = nxl_y_l, nxr_y_l
[1106]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
[1216]1143                      ar_tr(j,i,k) = work(2*j)
[1106]1144                   ENDDO
1145                   DO  j = 1, (ny+1)/2 - 1
[1216]1146                      ar_tr(ny+1-j,i,k) = work(2*j+1)
[1106]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
[1216]1158                DO  i = nxl_y_l, nxr_y_l
[1106]1159
1160                   DO  j = 0, (ny+1)/2
[1216]1161                      work(2*j) = ar_tr(j,i,k)
[1106]1162                   ENDDO
1163                   DO  j = 1, (ny+1)/2 - 1
[1216]1164                      work(2*j+1) = ar_tr(ny+1-j,i,k)
[1106]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
[1216]1187                DO  i = nxl_y_l, nxr_y_l
[1106]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
[1216]1194                      ar_tr(j,i,k) = work(2*j)
[1106]1195                   ENDDO
1196                   DO  j = 1, (ny+1)/2 - 1
[1216]1197                      ar_tr(ny+1-j,i,k) = work(2*j+1)
[1106]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
[1216]1209                DO  i = nxl_y_l, nxr_y_l
[1106]1210
1211                   DO  j = 0, (ny+1)/2
[1216]1212                      work(2*j) = ar_tr(j,i,k)
[1106]1213                   ENDDO
1214                   DO  j = 1, (ny+1)/2 - 1
[1216]1215                      work(2*j+1) = ar_tr(ny+1-j,i,k)
[1106]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
[1111]1233             !$acc data present( ar )
1234             CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
[1106]1235
[1111]1236             !$acc kernels
[1106]1237             DO  k = nzb_y, nzt_y
1238                DO  i = nxl_y, nxr_y
1239
1240                   DO  j = 0, (ny+1)/2
[1111]1241                      ar(j,i,k)      = REAL( ar_tmp(j,i,k) )  * dny
[1106]1242                   ENDDO
1243
1244                   DO  j = 1, (ny+1)/2 - 1
[1111]1245                      ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
[1106]1246                   ENDDO
1247
1248                ENDDO
1249             ENDDO
[1111]1250             !$acc end kernels
1251             !$acc end data
[1106]1252
1253          ELSE
1254
[1111]1255             !$acc data present( ar )
1256             !$acc kernels
[1106]1257             DO  k = nzb_y, nzt_y
1258                DO  i = nxl_y, nxr_y
1259
[1111]1260                   ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
[1106]1261
1262                   DO  j = 1, (ny+1)/2 - 1
[1111]1263                      ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
[1106]1264                   ENDDO
[1111]1265                   ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
[1106]1266
1267                ENDDO
1268             ENDDO
[1111]1269             !$acc end kernels
[1106]1270
[1111]1271             CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
1272             !$acc end data
[1106]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
[1]1306       INTEGER ::  j, jshape(1)
1307
[1106]1308       LOGICAL ::  forward_fft
1309
[1]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
[1106]1320       IF ( direction == 'forward' )  THEN
1321          forward_fft = .TRUE.
1322       ELSE
1323          forward_fft = .FALSE.
1324       ENDIF
1325
[1]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
[1106]1333          IF ( forward_fft )  THEN
[1]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
[1106]1374          IF ( forward_fft )  THEN
[1]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
[1216]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
[1]1431       ELSEIF ( fft_method == 'system-specific' )  THEN
1432
1433#if defined( __ibm )  &&  ! defined( __ibmy_special )
[1106]1434          IF ( forward_fft )  THEN
[1]1435
[1106]1436             CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
[1]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
[1106]1457             CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
[1]1458                         auy4, nau2 )
1459
1460             DO  j = 0, ny
1461                ar(j) = work(j)
1462             ENDDO
1463
1464          ENDIF
1465#elif defined( __nec )
[1106]1466          IF ( forward_fft )  THEN
[1]1467
1468             work(0:ny) = ar(0:ny)
1469
[1106]1470             CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
[1]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
[1106]1490             CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
[1]1491
1492             ar(0:ny) = work(0:ny)
1493
1494          ENDIF
1495#else
[254]1496          message_string = 'no system-specific fft-call available'
[1106]1497          CALL message( 'fft_y_1d', 'PA0188', 1, 2, 0, 6, 0 ) 
[254]1498
[1]1499#endif
1500
1501       ELSE
1502
[274]1503          message_string = 'fft method "' // TRIM( fft_method) // &
1504                           '" not available'
[1106]1505          CALL message( 'fft_y_1d', 'PA0189', 1, 2, 0, 6, 0 )
[1]1506
1507       ENDIF
1508
[1106]1509    END SUBROUTINE fft_y_1d
[1]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
[1092]1527       INTEGER                        ::  i, k, siza
[1]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 )
[1092]1533       INTEGER                             ::  sizw
[1]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
[1106]1587             CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
[1]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
[1106]1595             CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw, &
[1]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
[1106]1612             CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
[1]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
[1106]1626             CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, &
[1]1627                          trig_xb, work1, 0 )
1628
1629             ar(0:nx,1:nz) = ai(0:nx,1:nz)
1630
1631          ENDIF
1632
1633#else
[254]1634          message_string = 'no system-specific fft-call available'
1635          CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 ) 
[1]1636#endif
1637
1638       ELSE
1639
[274]1640          message_string = 'fft method "' // TRIM( fft_method) // &
1641                           '" not available'
[254]1642          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
[1]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
[1092]1664       INTEGER                        ::  j, k, ny1, siza
[1]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 )
[1092]1670       INTEGER                             ::  sizw
[1]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
[1106]1724             CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
[1]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
[1106]1732             CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, &
[1]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
[1106]1749             CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
[1]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
[1106]1763             CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, &
[1]1764                          trig_yb, work1, 0 )
1765
1766             ar(0:ny,1:nz) = ai(0:ny,1:nz)
1767
1768          ENDIF
1769
1770#else
[254]1771          message_string = 'no system-specific fft-call available'
1772          CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 ) 
[1]1773#endif
1774
1775       ELSE
[254]1776         
[274]1777          message_string = 'fft method "' // TRIM( fft_method) // &
1778                           '" not available'
[254]1779          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
[1]1780
1781       ENDIF
1782
1783    END SUBROUTINE fft_y_m
1784
[1106]1785
[1]1786 END MODULE fft_xy
Note: See TracBrowser for help on using the repository browser.