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

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

bugfixes:
missing variables added to ONLY lists in USE statements (advec_s_bc, advec_s_pw, advec_s_up, advec_ws, buoyancy, diffusion_e, diffusion_s, fft_xy, flow_statistics, palm, prognostic_equations)
missing module kinds added (cuda_fft_interfaces)
dpk renamed dp (fft_xy)
missing dependency for check_open added (Makefile)
variables removed from acc-present-list (diffusion_e, diffusion_w, diffusivities, production_e, wall_fluxes)
syntax errors removed from openacc-branch (flow_statistics)
USE-statement for nopointer-case added (swap_timelevel)

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