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

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