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

Last change on this file since 4516 was 4370, checked in by raasch, 5 years ago

bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added

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