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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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