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

Last change on this file since 4178 was 4069, checked in by Giersch, 2 years ago

Bugfix for masked output, compiler warning removed, test case for wind turbine model revised

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