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

Last change on this file since 4366 was 4366, checked in by raasch, 3 years ago

code vectorization for NEC Aurora: vectorized version of Temperton FFT, vectorization of Newtor iteration for calculating the Obukhov length

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