source: palm/trunk/SOURCE/fft_xy.f90 @ 1613

Last change on this file since 1613 was 1601, checked in by raasch, 9 years ago

last commit documented

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