source: palm/trunk/SOURCE/fft_xy_mod.f90

Last change on this file was 4828, checked in by Giersch, 6 months ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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