source: palm/trunk/SOURCE/fft_xy_mod.f90 @ 3634

Last change on this file since 3634 was 3634, checked in by knoop, 3 years ago

OpenACC port for SPEC

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