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

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

last commit documented

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