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

Last change on this file since 3774 was 3655, checked in by knoop, 6 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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