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

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

New:
---

GPU porting of pres, swap_timelevel. Adjustments of openACC directives.
Further porting of poisfft, which now runs completely on GPU without any
host/device data transfer for serial an parallel runs (but parallel runs
require data transfer before and after the MPI transpositions).
GPU-porting of tridiagonal solver:
tridiagonal routines split into extermal subroutines (instead using CONTAINS),
no distinction between parallel/non-parallel in poisfft and tridia any more,
tridia routines moved to end of file because of probable bug in PGI compiler
(otherwise "invalid device function" is indicated during runtime).
(cuda_fft_interfaces, fft_xy, flow_statistics, init_3d_model, palm, poisfft, pres, prognostic_equations, swap_timelevel, time_integration, transpose)
output of accelerator board information. (header)

optimization of tridia routines: constant elements and coefficients of tri are
stored in seperate arrays ddzuw and tric, last dimension of tri reduced from 5 to 2,
(init_grid, init_3d_model, modules, palm, poisfft)

poisfft_init is now called internally from poisfft,
(Makefile, Makefile_check, init_pegrid, poisfft, poisfft_hybrid)

CPU-time per grid point and timestep is output to CPU_MEASURES file
(cpu_statistics, modules, time_integration)

Changed:


resorting from/to array work changed, work now has 4 dimensions instead of 1 (transpose)
array diss allocated only if required (init_3d_model)

pressure boundary condition "Neumann+inhomo" removed from the code
(check_parameters, header, poisfft, poisfft_hybrid, pres)

Errors:


bugfix: dependency added for cuda_fft_interfaces (Makefile)
bugfix: CUDA fft plans adjusted for domain decomposition (before they always
used total domain) (fft_xy)

  • Property svn:keywords set to Id
File size: 45.6 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-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22! further openACC statements added, CUDA branch completely runs on GPU
23! bugfix: CUDA fft plans adjusted for domain decomposition (before they always
24! used total domain)
25!
26! Former revisions:
27! -----------------
28! $Id: fft_xy.f90 1111 2013-03-08 23:54:10Z raasch $
29!
30! 1106 2013-03-04 05:31:38Z raasch
31! CUDA fft added
32! array_kind renamed precision_kind, 3D- instead of 1D-loops in fft_x and fft_y
33! old fft_x, fft_y become fft_x_1d, fft_y_1d and are used for 1D-decomposition
34!
35! 1092 2013-02-02 11:24:22Z raasch
36! variable sizw declared for NEC case only
37!
38! 1036 2012-10-22 13:43:42Z raasch
39! code put under GPL (PALM 3.9)
40!
41! 274 2009-03-26 15:11:21Z heinze
42! Output of messages replaced by message handling routine.
43!
44! Feb. 2007
45! RCS Log replace by Id keyword, revision history cleaned up
46!
47! Revision 1.4  2006/03/28 12:27:09  raasch
48! Stop when system-specific fft is selected on NEC. For unknown reasons this
49! causes a program abort during first allocation in init_grid.
50!
51! Revision 1.2  2004/04/30 11:44:27  raasch
52! Module renamed from fft_for_1d_decomp to fft_xy, 1d-routines renamed to
53! fft_x and fft_y,
54! function FFT replaced by subroutine FFTN due to problems with 64-bit
55! mode on ibm,
56! shape of array cwork is explicitly stored in ishape/jshape and handled
57! to routine FFTN instead of shape-function (due to compiler error on
58! decalpha),
59! non vectorized FFT for nec included
60!
61! Revision 1.1  2002/06/11 13:00:49  raasch
62! Initial revision
63!
64!
65! Description:
66! ------------
67! Fast Fourier transformation along x and y for 1d domain decomposition along x.
68! Original version: Klaus Ketelsen (May 2002)
69!------------------------------------------------------------------------------!
70
71    USE control_parameters
72    USE indices
73    USE precision_kind
74    USE singleton
75    USE temperton_fft
76    USE transpose_indices
77
78    IMPLICIT NONE
79
80    PRIVATE
81    PUBLIC fft_x, fft_x_1d, fft_y, fft_y_1d, fft_init, fft_x_m, fft_y_m
82
83    INTEGER, DIMENSION(:), ALLOCATABLE, SAVE ::  ifax_x, ifax_y
84
85    LOGICAL, SAVE                            ::  init_fft = .FALSE.
86
87    REAL, SAVE ::  dnx, dny, sqr_dnx, sqr_dny
88    REAL, DIMENSION(:), ALLOCATABLE, SAVE    ::  trigs_x, trigs_y
89
90#if defined( __ibm )
91    INTEGER, PARAMETER ::  nau1 = 20000, nau2 = 22000
92!
93!-- The following working arrays contain tables and have to be "save" and
94!-- shared in OpenMP sense
95    REAL, DIMENSION(nau1), SAVE ::  aux1, auy1, aux3, auy3
96#elif defined( __nec )
97    INTEGER, SAVE ::  nz1
98    REAL, DIMENSION(:), ALLOCATABLE, SAVE ::  trig_xb, trig_xf, trig_yb, &
99                                              trig_yf
100#elif defined( __cuda_fft )
101    INTEGER, SAVE ::  plan_xf, plan_xi, plan_yf, plan_yi, total_points_x_transpo, &
102                      total_points_y_transpo
103#endif
104
105!
106!-- Public interfaces
107    INTERFACE fft_init
108       MODULE PROCEDURE fft_init
109    END INTERFACE fft_init
110
111    INTERFACE fft_x
112       MODULE PROCEDURE fft_x
113    END INTERFACE fft_x
114
115    INTERFACE fft_x_1d
116       MODULE PROCEDURE fft_x_1d
117    END INTERFACE fft_x_1d
118
119    INTERFACE fft_y
120       MODULE PROCEDURE fft_y
121    END INTERFACE fft_y
122
123    INTERFACE fft_y_1d
124       MODULE PROCEDURE fft_y_1d
125    END INTERFACE fft_y_1d
126
127    INTERFACE fft_x_m
128       MODULE PROCEDURE fft_x_m
129    END INTERFACE fft_x_m
130
131    INTERFACE fft_y_m
132       MODULE PROCEDURE fft_y_m
133    END INTERFACE fft_y_m
134
135 CONTAINS
136
137
138    SUBROUTINE fft_init
139
140       USE cuda_fft_interfaces
141
142       IMPLICIT NONE
143
144!
145!--    The following temporary working arrays have to be on stack or private
146!--    in OpenMP sense
147#if defined( __ibm )
148       REAL, DIMENSION(0:nx+2) :: workx
149       REAL, DIMENSION(0:ny+2) :: worky
150       REAL, DIMENSION(nau2)   :: aux2, auy2, aux4, auy4
151#elif defined( __nec )
152       REAL, DIMENSION(0:nx+3,nz+1)   ::  work_x
153       REAL, DIMENSION(0:ny+3,nz+1)   ::  work_y
154       REAL, DIMENSION(6*(nx+3),nz+1) ::  workx
155       REAL, DIMENSION(6*(ny+3),nz+1) ::  worky
156#endif
157
158!
159!--    Return, if already called
160       IF ( init_fft )  THEN
161          RETURN
162       ELSE
163          init_fft = .TRUE.
164       ENDIF
165
166       IF ( fft_method == 'system-specific' )  THEN
167
168          dnx = 1.0 / ( nx + 1.0 )
169          dny = 1.0 / ( ny + 1.0 )
170          sqr_dnx = SQRT( dnx )
171          sqr_dny = SQRT( dny )
172#if defined( __ibm )  &&  ! defined( __ibmy_special )
173!
174!--       Initialize tables for fft along x
175          CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1,  1, sqr_dnx, aux1, nau1, &
176                      aux2, nau2 )
177          CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
178                      aux4, nau2 )
179!
180!--       Initialize tables for fft along y
181          CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1,  1, sqr_dny, auy1, nau1, &
182                      auy2, nau2 )
183          CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
184                      auy4, nau2 )
185#elif defined( __nec )
186          message_string = 'fft method "' // TRIM( fft_method) // &
187                           '" currently does not work on NEC'
188          CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 )
189
190          ALLOCATE( trig_xb(2*(nx+1)), trig_xf(2*(nx+1)), &
191                    trig_yb(2*(ny+1)), trig_yf(2*(ny+1)) )
192
193          work_x = 0.0
194          work_y = 0.0
195          nz1  = nz + MOD( nz+1, 2 )  ! odd nz slows down fft significantly
196                                      ! when using the NEC ffts
197
198!
199!--       Initialize tables for fft along x (non-vector and vector case (M))
200          CALL DZFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xf, workx, 0 )
201          CALL ZDFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xb, workx, 0 )
202          CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
203                       trig_xf, workx, 0 )
204          CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
205                       trig_xb, workx, 0 )
206!
207!--       Initialize tables for fft along y (non-vector and vector case (M))
208          CALL DZFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yf, worky, 0 )
209          CALL ZDFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yb, worky, 0 )
210          CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
211                       trig_yf, worky, 0 )
212          CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
213                       trig_yb, worky, 0 )
214#elif defined( __cuda_fft )
215          total_points_x_transpo = (nx+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1)
216          total_points_y_transpo = (ny+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1)
217          CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
218          CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
219          CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
220          CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
221#else
222          message_string = 'no system-specific fft-call available'
223          CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
224#endif
225       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
226!
227!--       Temperton-algorithm
228!--       Initialize tables for fft along x and y
229          ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
230
231          CALL set99( trigs_x, ifax_x, nx+1 )
232          CALL set99( trigs_y, ifax_y, ny+1 )
233
234       ELSEIF ( fft_method == 'singleton-algorithm' )  THEN
235
236          CONTINUE
237
238       ELSE
239
240          message_string = 'fft method "' // TRIM( fft_method) // &
241                           '" not available'
242          CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
243       ENDIF
244
245    END SUBROUTINE fft_init
246
247
248    SUBROUTINE fft_x( ar, direction )
249
250!----------------------------------------------------------------------!
251!                                 fft_x                                !
252!                                                                      !
253!               Fourier-transformation along x-direction               !
254!                     Version for 2D-decomposition                     !
255!                                                                      !
256!      fft_x uses internal algorithms (Singleton or Temperton) or      !
257!           system-specific routines, if they are available            !
258!----------------------------------------------------------------------!
259
260       USE cuda_fft_interfaces
261
262       IMPLICIT NONE
263
264       CHARACTER (LEN=*) ::  direction
265       INTEGER ::  i, ishape(1), j, k
266
267       LOGICAL ::  forward_fft
268
269       REAL, DIMENSION(0:nx+2)   ::  work
270       REAL, DIMENSION(nx+2)     ::  work1
271       COMPLEX, DIMENSION(:), ALLOCATABLE ::  cwork
272#if defined( __ibm )
273       REAL, DIMENSION(nau2)     ::  aux2, aux4
274#elif defined( __nec )
275       REAL, DIMENSION(6*(nx+1)) ::  work2
276#elif defined( __cuda_fft )
277       !$acc declare create( ar_tmp )
278       COMPLEX(dpk), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) ::  ar_tmp
279#endif
280       REAL, DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::  ar
281
282       IF ( direction == 'forward' )  THEN
283          forward_fft = .TRUE.
284       ELSE
285          forward_fft = .FALSE.
286       ENDIF
287
288       IF ( fft_method == 'singleton-algorithm' )  THEN
289
290!
291!--       Performing the fft with singleton's software works on every system,
292!--       since it is part of the model
293          ALLOCATE( cwork(0:nx) )
294     
295          IF ( forward_fft )   then
296
297             !$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
298             !$OMP DO
299             DO  k = nzb_x, nzt_x
300                DO  j = nys_x, nyn_x
301
302                   DO  i = 0, nx
303                      cwork(i) = CMPLX( ar(i,j,k) )
304                   ENDDO
305
306                   ishape = SHAPE( cwork )
307                   CALL FFTN( cwork, ishape )
308
309                   DO  i = 0, (nx+1)/2
310                      ar(i,j,k) = REAL( cwork(i) )
311                   ENDDO
312                   DO  i = 1, (nx+1)/2 - 1
313                      ar(nx+1-i,j,k) = -AIMAG( cwork(i) )
314                   ENDDO
315
316                ENDDO
317             ENDDO
318             !$OMP END PARALLEL
319
320          ELSE
321
322             !$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
323             !$OMP DO
324             DO  k = nzb_x, nzt_x
325                DO  j = nys_x, nyn_x
326
327                   cwork(0) = CMPLX( ar(0,j,k), 0.0 )
328                   DO  i = 1, (nx+1)/2 - 1
329                      cwork(i)      = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k) )
330                      cwork(nx+1-i) = CMPLX( ar(i,j,k),  ar(nx+1-i,j,k) )
331                   ENDDO
332                   cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
333
334                   ishape = SHAPE( cwork )
335                   CALL FFTN( cwork, ishape, inv = .TRUE. )
336
337                   DO  i = 0, nx
338                      ar(i,j,k) = REAL( cwork(i) )
339                   ENDDO
340
341                ENDDO
342             ENDDO
343             !$OMP END PARALLEL
344
345          ENDIF
346
347          DEALLOCATE( cwork )
348
349       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
350
351!
352!--       Performing the fft with Temperton's software works on every system,
353!--       since it is part of the model
354          IF ( forward_fft )  THEN
355
356             !$OMP PARALLEL PRIVATE ( work, i, j, k )
357             !$OMP DO
358             DO  k = nzb_x, nzt_x
359                DO  j = nys_x, nyn_x
360
361                   work(0:nx) = ar(0:nx,j,k)
362                   CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
363
364                   DO  i = 0, (nx+1)/2
365                      ar(i,j,k) = work(2*i)
366                   ENDDO
367                   DO  i = 1, (nx+1)/2 - 1
368                      ar(nx+1-i,j,k) = work(2*i+1)
369                   ENDDO
370
371                ENDDO
372             ENDDO
373             !$OMP END PARALLEL
374
375          ELSE
376
377             !$OMP PARALLEL PRIVATE ( work, i, j, k )
378             !$OMP DO
379             DO  k = nzb_x, nzt_x
380                DO  j = nys_x, nyn_x
381
382                   DO  i = 0, (nx+1)/2
383                      work(2*i) = ar(i,j,k)
384                   ENDDO
385                   DO  i = 1, (nx+1)/2 - 1
386                      work(2*i+1) = ar(nx+1-i,j,k)
387                   ENDDO
388                   work(1)    = 0.0
389                   work(nx+2) = 0.0
390
391                   CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
392                   ar(0:nx,j,k) = work(0:nx)
393
394                ENDDO
395             ENDDO
396             !$OMP END PARALLEL
397
398          ENDIF
399
400       ELSEIF ( fft_method == 'system-specific' )  THEN
401
402#if defined( __ibm )  &&  ! defined( __ibmy_special )
403          IF ( forward_fft )  THEN
404
405             !$OMP PARALLEL PRIVATE ( work, i, j, k )
406             !$OMP DO
407             DO  k = nzb_x, nzt_x
408                DO  j = nys_x, nyn_x
409
410                   CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
411                               aux2, nau2 )
412
413                   DO  i = 0, (nx+1)/2
414                      ar(i,j,k) = work(2*i)
415                   ENDDO
416                   DO  i = 1, (nx+1)/2 - 1
417                      ar(nx+1-i,j,k) = work(2*i+1)
418                   ENDDO
419
420                ENDDO
421             ENDDO
422             !$OMP END PARALLEL
423
424          ELSE
425
426             !$OMP PARALLEL PRIVATE ( work, i, j, k )
427             !$OMP DO
428             DO  k = nzb_x, nzt_x
429                DO  j = nys_x, nyn_x
430
431                   DO  i = 0, (nx+1)/2
432                      work(2*i) = ar(i,j,k)
433                   ENDDO
434                   DO  i = 1, (nx+1)/2 - 1
435                      work(2*i+1) = ar(nx+1-i,j,k)
436                   ENDDO
437                   work(1) = 0.0
438                   work(nx+2) = 0.0
439
440                   CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
441                               aux4, nau2 )
442
443                   DO  i = 0, nx
444                      ar(i,j,k) = work(i)
445                   ENDDO
446
447                ENDDO
448             ENDDO
449             !$OMP END PARALLEL
450
451          ENDIF
452
453#elif defined( __nec )
454
455          IF ( forward_fft )  THEN
456
457             !$OMP PARALLEL PRIVATE ( work, i, j, k )
458             !$OMP DO
459             DO  k = nzb_x, nzt_x
460                DO  j = nys_x, nyn_x
461
462                   work(0:nx) = ar(0:nx,j,k)
463
464                   CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
465     
466                   DO  i = 0, (nx+1)/2
467                      ar(i,j,k) = work(2*i)
468                   ENDDO
469                   DO  i = 1, (nx+1)/2 - 1
470                      ar(nx+1-i,j,k) = work(2*i+1)
471                   ENDDO
472
473                ENDDO
474             ENDDO
475             !$END OMP PARALLEL
476
477          ELSE
478
479             !$OMP PARALLEL PRIVATE ( work, i, j, k )
480             !$OMP DO
481             DO  k = nzb_x, nzt_x
482                DO  j = nys_x, nyn_x
483
484                   DO  i = 0, (nx+1)/2
485                      work(2*i) = ar(i,j,k)
486                   ENDDO
487                   DO  i = 1, (nx+1)/2 - 1
488                      work(2*i+1) = ar(nx+1-i,j,k)
489                   ENDDO
490                   work(1) = 0.0
491                   work(nx+2) = 0.0
492
493                   CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
494
495                   ar(0:nx,j,k) = work(0:nx)
496
497                ENDDO
498             ENDDO
499             !$OMP END PARALLEL
500
501          ENDIF
502
503#elif defined( __cuda_fft )
504
505          IF ( forward_fft )  THEN
506
507             !$acc data present( ar )
508             CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
509
510             !$acc kernels
511             !$acc loop
512             DO  k = nzb_x, nzt_x
513                DO  j = nys_x, nyn_x
514
515                   !$acc loop vector( 32 )
516                   DO  i = 0, (nx+1)/2
517                      ar(i,j,k)      = REAL( ar_tmp(i,j,k) )  * dnx
518                   ENDDO
519
520                   !$acc loop vector( 32 )
521                   DO  i = 1, (nx+1)/2 - 1
522                      ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
523                   ENDDO
524
525                ENDDO
526             ENDDO
527             !$acc end kernels
528             !$acc end data
529
530          ELSE
531
532             !$acc data present( ar )
533             !$acc kernels
534             !$acc loop
535             DO  k = nzb_x, nzt_x
536                DO  j = nys_x, nyn_x
537
538                   ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
539
540                   !$acc loop vector( 32 )
541                   DO  i = 1, (nx+1)/2 - 1
542                      ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
543                   ENDDO
544                   ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
545
546                ENDDO
547             ENDDO
548             !$acc end kernels
549
550             CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
551             !$acc end data
552
553          ENDIF
554
555#else
556          message_string = 'no system-specific fft-call available'
557          CALL message( 'fft_x', 'PA0188', 1, 2, 0, 6, 0 )
558#endif
559
560       ELSE
561
562          message_string = 'fft method "' // TRIM( fft_method) // &
563                           '" not available'
564          CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
565
566       ENDIF
567
568    END SUBROUTINE fft_x
569
570    SUBROUTINE fft_x_1d( ar, direction )
571
572!----------------------------------------------------------------------!
573!                               fft_x_1d                               !
574!                                                                      !
575!               Fourier-transformation along x-direction               !
576!                     Version for 1D-decomposition                     !
577!                                                                      !
578!      fft_x uses internal algorithms (Singleton or Temperton) or      !
579!           system-specific routines, if they are available            !
580!----------------------------------------------------------------------!
581
582       IMPLICIT NONE
583
584       CHARACTER (LEN=*) ::  direction
585       INTEGER ::  i, ishape(1)
586
587       LOGICAL ::  forward_fft
588
589       REAL, DIMENSION(0:nx)     ::  ar
590       REAL, DIMENSION(0:nx+2)   ::  work
591       REAL, DIMENSION(nx+2)     ::  work1
592       COMPLEX, DIMENSION(:), ALLOCATABLE ::  cwork
593#if defined( __ibm )
594       REAL, DIMENSION(nau2)     ::  aux2, aux4
595#elif defined( __nec )
596       REAL, DIMENSION(6*(nx+1)) ::  work2
597#endif
598
599       IF ( direction == 'forward' )  THEN
600          forward_fft = .TRUE.
601       ELSE
602          forward_fft = .FALSE.
603       ENDIF
604
605       IF ( fft_method == 'singleton-algorithm' )  THEN
606
607!
608!--       Performing the fft with singleton's software works on every system,
609!--       since it is part of the model
610          ALLOCATE( cwork(0:nx) )
611     
612          IF ( forward_fft )   then
613
614             DO  i = 0, nx
615                cwork(i) = CMPLX( ar(i) )
616             ENDDO
617             ishape = SHAPE( cwork )
618             CALL FFTN( cwork, ishape )
619             DO  i = 0, (nx+1)/2
620                ar(i) = REAL( cwork(i) )
621             ENDDO
622             DO  i = 1, (nx+1)/2 - 1
623                ar(nx+1-i) = -AIMAG( cwork(i) )
624             ENDDO
625
626          ELSE
627
628             cwork(0) = CMPLX( ar(0), 0.0 )
629             DO  i = 1, (nx+1)/2 - 1
630                cwork(i)      = CMPLX( ar(i), -ar(nx+1-i) )
631                cwork(nx+1-i) = CMPLX( ar(i),  ar(nx+1-i) )
632             ENDDO
633             cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
634
635             ishape = SHAPE( cwork )
636             CALL FFTN( cwork, ishape, inv = .TRUE. )
637
638             DO  i = 0, nx
639                ar(i) = REAL( cwork(i) )
640             ENDDO
641
642          ENDIF
643
644          DEALLOCATE( cwork )
645
646       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
647
648!
649!--       Performing the fft with Temperton's software works on every system,
650!--       since it is part of the model
651          IF ( forward_fft )  THEN
652
653             work(0:nx) = ar
654             CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
655
656             DO  i = 0, (nx+1)/2
657                ar(i) = work(2*i)
658             ENDDO
659             DO  i = 1, (nx+1)/2 - 1
660                ar(nx+1-i) = work(2*i+1)
661             ENDDO
662
663          ELSE
664
665             DO  i = 0, (nx+1)/2
666                work(2*i) = ar(i)
667             ENDDO
668             DO  i = 1, (nx+1)/2 - 1
669                work(2*i+1) = ar(nx+1-i)
670             ENDDO
671             work(1)    = 0.0
672             work(nx+2) = 0.0
673
674             CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
675             ar = work(0:nx)
676
677          ENDIF
678
679       ELSEIF ( fft_method == 'system-specific' )  THEN
680
681#if defined( __ibm )  &&  ! defined( __ibmy_special )
682          IF ( forward_fft )  THEN
683
684             CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
685                         aux2, nau2 )
686
687             DO  i = 0, (nx+1)/2
688                ar(i) = work(2*i)
689             ENDDO
690             DO  i = 1, (nx+1)/2 - 1
691                ar(nx+1-i) = work(2*i+1)
692             ENDDO
693
694          ELSE
695
696             DO  i = 0, (nx+1)/2
697                work(2*i) = ar(i)
698             ENDDO
699             DO  i = 1, (nx+1)/2 - 1
700                work(2*i+1) = ar(nx+1-i)
701             ENDDO
702             work(1) = 0.0
703             work(nx+2) = 0.0
704
705             CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
706                         aux4, nau2 )
707
708             DO  i = 0, nx
709                ar(i) = work(i)
710             ENDDO
711
712          ENDIF
713#elif defined( __nec )
714          IF ( forward_fft )  THEN
715
716             work(0:nx) = ar(0:nx)
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) = work(2*i)
722             ENDDO
723             DO  i = 1, (nx+1)/2 - 1
724                ar(nx+1-i) = work(2*i+1)
725             ENDDO
726
727          ELSE
728
729             DO  i = 0, (nx+1)/2
730                work(2*i) = ar(i)
731             ENDDO
732             DO  i = 1, (nx+1)/2 - 1
733                work(2*i+1) = ar(nx+1-i)
734             ENDDO
735             work(1) = 0.0
736             work(nx+2) = 0.0
737
738             CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
739
740             ar(0:nx) = work(0:nx)
741
742          ENDIF
743#else
744          message_string = 'no system-specific fft-call available'
745          CALL message( 'fft_x_1d', 'PA0188', 1, 2, 0, 6, 0 )
746#endif
747       ELSE
748          message_string = 'fft method "' // TRIM( fft_method) // &
749                           '" not available'
750          CALL message( 'fft_x_1d', 'PA0189', 1, 2, 0, 6, 0 )
751
752       ENDIF
753
754    END SUBROUTINE fft_x_1d
755
756    SUBROUTINE fft_y( ar, direction )
757
758!----------------------------------------------------------------------!
759!                                 fft_y                                !
760!                                                                      !
761!               Fourier-transformation along y-direction               !
762!                     Version for 2D-decomposition                     !
763!                                                                      !
764!      fft_y uses internal algorithms (Singleton or Temperton) or      !
765!           system-specific routines, if they are available            !
766!----------------------------------------------------------------------!
767
768       USE cuda_fft_interfaces
769
770       IMPLICIT NONE
771
772       CHARACTER (LEN=*) ::  direction
773       INTEGER ::  i, j, jshape(1), k
774
775       LOGICAL ::  forward_fft
776
777       REAL, DIMENSION(0:ny+2)   ::  work
778       REAL, DIMENSION(ny+2)     ::  work1
779       COMPLEX, DIMENSION(:), ALLOCATABLE ::  cwork
780#if defined( __ibm )
781       REAL, DIMENSION(nau2)     ::  auy2, auy4
782#elif defined( __nec )
783       REAL, DIMENSION(6*(ny+1)) ::  work2
784#elif defined( __cuda_fft )
785       !$acc declare create( ar_tmp )
786       COMPLEX(dpk), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) ::  ar_tmp
787#endif
788       REAL, DIMENSION(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) ::  ar
789
790       IF ( direction == 'forward' )  THEN
791          forward_fft = .TRUE.
792       ELSE
793          forward_fft = .FALSE.
794       ENDIF
795
796       IF ( fft_method == 'singleton-algorithm' )  THEN
797
798!
799!--       Performing the fft with singleton's software works on every system,
800!--       since it is part of the model
801          ALLOCATE( cwork(0:ny) )
802
803          IF ( forward_fft )   then
804
805             !$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
806             !$OMP DO
807             DO  k = nzb_y, nzt_y
808                DO  i = nxl_y, nxr_y
809
810                   DO  j = 0, ny
811                      cwork(j) = CMPLX( ar(j,i,k) )
812                   ENDDO
813
814                   jshape = SHAPE( cwork )
815                   CALL FFTN( cwork, jshape )
816
817                   DO  j = 0, (ny+1)/2
818                      ar(j,i,k) = REAL( cwork(j) )
819                   ENDDO
820                   DO  j = 1, (ny+1)/2 - 1
821                      ar(ny+1-j,i,k) = -AIMAG( cwork(j) )
822                   ENDDO
823
824                ENDDO
825             ENDDO
826             !$OMP END PARALLEL
827
828          ELSE
829
830             !$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
831             !$OMP DO
832             DO  k = nzb_y, nzt_y
833                DO  i = nxl_y, nxr_y
834
835                   cwork(0) = CMPLX( ar(0,i,k), 0.0 )
836                   DO  j = 1, (ny+1)/2 - 1
837                      cwork(j)      = CMPLX( ar(j,i,k), -ar(ny+1-j,i,k) )
838                      cwork(ny+1-j) = CMPLX( ar(j,i,k),  ar(ny+1-j,i,k) )
839                   ENDDO
840                   cwork((ny+1)/2) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
841
842                   jshape = SHAPE( cwork )
843                   CALL FFTN( cwork, jshape, inv = .TRUE. )
844
845                   DO  j = 0, ny
846                      ar(j,i,k) = REAL( cwork(j) )
847                   ENDDO
848
849                ENDDO
850             ENDDO
851             !$OMP END PARALLEL
852
853          ENDIF
854
855          DEALLOCATE( cwork )
856
857       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
858
859!
860!--       Performing the fft with Temperton's software works on every system,
861!--       since it is part of the model
862          IF ( forward_fft )  THEN
863
864             !$OMP PARALLEL PRIVATE ( work, i, j, k )
865             !$OMP DO
866             DO  k = nzb_y, nzt_y
867                DO  i = nxl_y, nxr_y
868
869                   work(0:ny) = ar(0:ny,i,k)
870                   CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
871
872                   DO  j = 0, (ny+1)/2
873                      ar(j,i,k) = work(2*j)
874                   ENDDO
875                   DO  j = 1, (ny+1)/2 - 1
876                      ar(ny+1-j,i,k) = work(2*j+1)
877                   ENDDO
878
879                ENDDO
880             ENDDO
881             !$OMP END PARALLEL
882
883          ELSE
884
885             !$OMP PARALLEL PRIVATE ( work, i, j, k )
886             !$OMP DO
887             DO  k = nzb_y, nzt_y
888                DO  i = nxl_y, nxr_y
889
890                   DO  j = 0, (ny+1)/2
891                      work(2*j) = ar(j,i,k)
892                   ENDDO
893                   DO  j = 1, (ny+1)/2 - 1
894                      work(2*j+1) = ar(ny+1-j,i,k)
895                   ENDDO
896                   work(1)    = 0.0
897                   work(ny+2) = 0.0
898
899                   CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
900                   ar(0:ny,i,k) = work(0:ny)
901
902                ENDDO
903             ENDDO
904             !$OMP END PARALLEL
905
906          ENDIF
907
908       ELSEIF ( fft_method == 'system-specific' )  THEN
909
910#if defined( __ibm )  &&  ! defined( __ibmy_special )
911          IF ( forward_fft)  THEN
912
913             !$OMP PARALLEL PRIVATE ( work, i, j, k )
914             !$OMP DO
915             DO  k = nzb_y, nzt_y
916                DO  i = nxl_y, nxr_y
917
918                   CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
919                               auy2, nau2 )
920
921                   DO  j = 0, (ny+1)/2
922                      ar(j,i,k) = work(2*j)
923                   ENDDO
924                   DO  j = 1, (ny+1)/2 - 1
925                      ar(ny+1-j,i,k) = work(2*j+1)
926                   ENDDO
927
928                ENDDO
929             ENDDO
930             !$OMP END PARALLEL
931
932          ELSE
933
934             !$OMP PARALLEL PRIVATE ( work, i, j, k )
935             !$OMP DO
936             DO  k = nzb_y, nzt_y
937                DO  i = nxl_y, nxr_y
938
939                   DO  j = 0, (ny+1)/2
940                      work(2*j) = ar(j,i,k)
941                   ENDDO
942                   DO  j = 1, (ny+1)/2 - 1
943                      work(2*j+1) = ar(ny+1-j,i,k)
944                   ENDDO
945                   work(1)    = 0.0
946                   work(ny+2) = 0.0
947
948                   CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
949                               auy4, nau2 )
950
951                   DO  j = 0, ny
952                      ar(j,i,k) = work(j)
953                   ENDDO
954
955                ENDDO
956             ENDDO
957             !$OMP END PARALLEL
958
959          ENDIF
960#elif defined( __nec )
961          IF ( forward_fft )  THEN
962
963             !$OMP PARALLEL PRIVATE ( work, i, j, k )
964             !$OMP DO
965             DO  k = nzb_y, nzt_y
966                DO  i = nxl_y, nxr_y
967
968                   work(0:ny) = ar(0:ny,i,k)
969
970                   CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
971
972                   DO  j = 0, (ny+1)/2
973                      ar(j,i,k) = work(2*j)
974                   ENDDO
975                   DO  j = 1, (ny+1)/2 - 1
976                      ar(ny+1-j,i,k) = work(2*j+1)
977                   ENDDO
978
979                ENDDO
980             ENDDO
981             !$END OMP PARALLEL
982
983          ELSE
984
985             !$OMP PARALLEL PRIVATE ( work, i, j, k )
986             !$OMP DO
987             DO  k = nzb_y, nzt_y
988                DO  i = nxl_y, nxr_y
989
990                   DO  j = 0, (ny+1)/2
991                      work(2*j) = ar(j,i,k)
992                   ENDDO
993                   DO  j = 1, (ny+1)/2 - 1
994                      work(2*j+1) = ar(ny+1-j,i,k)
995                   ENDDO
996                   work(1) = 0.0
997                   work(ny+2) = 0.0
998
999                   CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1000
1001                   ar(0:ny,i,k) = work(0:ny)
1002
1003                ENDDO
1004             ENDDO
1005             !$OMP END PARALLEL
1006
1007          ENDIF
1008#elif defined( __cuda_fft )
1009
1010          IF ( forward_fft )  THEN
1011
1012             !$acc data present( ar )
1013             CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
1014
1015             !$acc kernels
1016             !$acc loop
1017             DO  k = nzb_y, nzt_y
1018                DO  i = nxl_y, nxr_y
1019
1020                   !$acc loop vector( 32 )
1021                   DO  j = 0, (ny+1)/2
1022                      ar(j,i,k)      = REAL( ar_tmp(j,i,k) )  * dny
1023                   ENDDO
1024
1025                   !$acc loop vector( 32 )
1026                   DO  j = 1, (ny+1)/2 - 1
1027                      ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
1028                   ENDDO
1029
1030                ENDDO
1031             ENDDO
1032             !$acc end kernels
1033             !$acc end data
1034
1035          ELSE
1036
1037             !$acc data present( ar )
1038             !$acc kernels
1039             !$acc loop
1040             DO  k = nzb_y, nzt_y
1041                DO  i = nxl_y, nxr_y
1042
1043                   ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
1044
1045                   !$acc loop vector( 32 )
1046                   DO  j = 1, (ny+1)/2 - 1
1047                      ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
1048                   ENDDO
1049                   ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
1050
1051                ENDDO
1052             ENDDO
1053             !$acc end kernels
1054
1055             CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
1056             !$acc end data
1057
1058          ENDIF
1059
1060#else
1061          message_string = 'no system-specific fft-call available'
1062          CALL message( 'fft_y', 'PA0188', 1, 2, 0, 6, 0 ) 
1063#endif
1064
1065       ELSE
1066
1067          message_string = 'fft method "' // TRIM( fft_method) // &
1068                           '" not available'
1069          CALL message( 'fft_y', 'PA0189', 1, 2, 0, 6, 0 )
1070
1071       ENDIF
1072
1073    END SUBROUTINE fft_y
1074
1075    SUBROUTINE fft_y_1d( ar, direction )
1076
1077!----------------------------------------------------------------------!
1078!                               fft_y_1d                               !
1079!                                                                      !
1080!               Fourier-transformation along y-direction               !
1081!                     Version for 1D-decomposition                     !
1082!                                                                      !
1083!      fft_y uses internal algorithms (Singleton or Temperton) or      !
1084!           system-specific routines, if they are available            !
1085!----------------------------------------------------------------------!
1086
1087       IMPLICIT NONE
1088
1089       CHARACTER (LEN=*) ::  direction
1090       INTEGER ::  j, jshape(1)
1091
1092       LOGICAL ::  forward_fft
1093
1094       REAL, DIMENSION(0:ny)     ::  ar
1095       REAL, DIMENSION(0:ny+2)   ::  work
1096       REAL, DIMENSION(ny+2)     ::  work1
1097       COMPLEX, DIMENSION(:), ALLOCATABLE ::  cwork
1098#if defined( __ibm )
1099       REAL, DIMENSION(nau2)     ::  auy2, auy4
1100#elif defined( __nec )
1101       REAL, DIMENSION(6*(ny+1)) ::  work2
1102#endif
1103
1104       IF ( direction == 'forward' )  THEN
1105          forward_fft = .TRUE.
1106       ELSE
1107          forward_fft = .FALSE.
1108       ENDIF
1109
1110       IF ( fft_method == 'singleton-algorithm' )  THEN
1111
1112!
1113!--       Performing the fft with singleton's software works on every system,
1114!--       since it is part of the model
1115          ALLOCATE( cwork(0:ny) )
1116
1117          IF ( forward_fft )  THEN
1118
1119             DO  j = 0, ny
1120                cwork(j) = CMPLX( ar(j) )
1121             ENDDO
1122
1123             jshape = SHAPE( cwork )
1124             CALL FFTN( cwork, jshape )
1125
1126             DO  j = 0, (ny+1)/2
1127                ar(j) = REAL( cwork(j) )
1128             ENDDO
1129             DO  j = 1, (ny+1)/2 - 1
1130                ar(ny+1-j) = -AIMAG( cwork(j) )
1131             ENDDO
1132
1133          ELSE
1134
1135             cwork(0) = CMPLX( ar(0), 0.0 )
1136             DO  j = 1, (ny+1)/2 - 1
1137                cwork(j)      = CMPLX( ar(j), -ar(ny+1-j) )
1138                cwork(ny+1-j) = CMPLX( ar(j),  ar(ny+1-j) )
1139             ENDDO
1140             cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
1141
1142             jshape = SHAPE( cwork )
1143             CALL FFTN( cwork, jshape, inv = .TRUE. )
1144
1145             DO  j = 0, ny
1146                ar(j) = REAL( cwork(j) )
1147             ENDDO
1148
1149          ENDIF
1150
1151          DEALLOCATE( cwork )
1152
1153       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
1154
1155!
1156!--       Performing the fft with Temperton's software works on every system,
1157!--       since it is part of the model
1158          IF ( forward_fft )  THEN
1159
1160             work(0:ny) = ar
1161             CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
1162
1163             DO  j = 0, (ny+1)/2
1164                ar(j) = work(2*j)
1165             ENDDO
1166             DO  j = 1, (ny+1)/2 - 1
1167                ar(ny+1-j) = work(2*j+1)
1168             ENDDO
1169
1170          ELSE
1171
1172             DO  j = 0, (ny+1)/2
1173                work(2*j) = ar(j)
1174             ENDDO
1175             DO  j = 1, (ny+1)/2 - 1
1176                work(2*j+1) = ar(ny+1-j)
1177             ENDDO
1178             work(1)    = 0.0
1179             work(ny+2) = 0.0
1180
1181             CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
1182             ar = work(0:ny)
1183
1184          ENDIF
1185
1186       ELSEIF ( fft_method == 'system-specific' )  THEN
1187
1188#if defined( __ibm )  &&  ! defined( __ibmy_special )
1189          IF ( forward_fft )  THEN
1190
1191             CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
1192                         auy2, nau2 )
1193
1194             DO  j = 0, (ny+1)/2
1195                ar(j) = work(2*j)
1196             ENDDO
1197             DO  j = 1, (ny+1)/2 - 1
1198                ar(ny+1-j) = work(2*j+1)
1199             ENDDO
1200
1201          ELSE
1202
1203             DO  j = 0, (ny+1)/2
1204                work(2*j) = ar(j)
1205             ENDDO
1206             DO  j = 1, (ny+1)/2 - 1
1207                work(2*j+1) = ar(ny+1-j)
1208             ENDDO
1209             work(1)    = 0.0
1210             work(ny+2) = 0.0
1211
1212             CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
1213                         auy4, nau2 )
1214
1215             DO  j = 0, ny
1216                ar(j) = work(j)
1217             ENDDO
1218
1219          ENDIF
1220#elif defined( __nec )
1221          IF ( forward_fft )  THEN
1222
1223             work(0:ny) = ar(0:ny)
1224
1225             CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
1226
1227             DO  j = 0, (ny+1)/2
1228                ar(j) = work(2*j)
1229             ENDDO
1230             DO  j = 1, (ny+1)/2 - 1
1231                ar(ny+1-j) = work(2*j+1)
1232             ENDDO
1233
1234          ELSE
1235
1236             DO  j = 0, (ny+1)/2
1237                work(2*j) = ar(j)
1238             ENDDO
1239             DO  j = 1, (ny+1)/2 - 1
1240                work(2*j+1) = ar(ny+1-j)
1241             ENDDO
1242             work(1) = 0.0
1243             work(ny+2) = 0.0
1244
1245             CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1246
1247             ar(0:ny) = work(0:ny)
1248
1249          ENDIF
1250#else
1251          message_string = 'no system-specific fft-call available'
1252          CALL message( 'fft_y_1d', 'PA0188', 1, 2, 0, 6, 0 ) 
1253
1254#endif
1255
1256       ELSE
1257
1258          message_string = 'fft method "' // TRIM( fft_method) // &
1259                           '" not available'
1260          CALL message( 'fft_y_1d', 'PA0189', 1, 2, 0, 6, 0 )
1261
1262       ENDIF
1263
1264    END SUBROUTINE fft_y_1d
1265
1266    SUBROUTINE fft_x_m( ar, direction )
1267
1268!----------------------------------------------------------------------!
1269!                               fft_x_m                                !
1270!                                                                      !
1271!               Fourier-transformation along x-direction               !
1272!                 Version for 1d domain decomposition                  !
1273!             using multiple 1D FFT from Math Keisan on NEC            !
1274!                       or Temperton-algorithm                         !
1275!       (no singleton-algorithm on NEC because it does not vectorize)  !
1276!                                                                      !
1277!----------------------------------------------------------------------!
1278
1279       IMPLICIT NONE
1280
1281       CHARACTER (LEN=*)              ::  direction
1282       INTEGER                        ::  i, k, siza
1283
1284       REAL, DIMENSION(0:nx,nz)       ::  ar
1285       REAL, DIMENSION(0:nx+3,nz+1)   ::  ai
1286       REAL, DIMENSION(6*(nx+4),nz+1) ::  work1
1287#if defined( __nec )
1288       INTEGER                             ::  sizw
1289       COMPLEX, DIMENSION((nx+4)/2+1,nz+1) ::  work
1290#endif
1291
1292       IF ( fft_method == 'temperton-algorithm' )  THEN
1293
1294          siza = SIZE( ai, 1 )
1295
1296          IF ( direction == 'forward')  THEN
1297
1298             ai(0:nx,1:nz) = ar(0:nx,1:nz)
1299             ai(nx+1:,:)   = 0.0
1300
1301             CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
1302
1303             DO  k = 1, nz
1304                DO  i = 0, (nx+1)/2
1305                   ar(i,k) = ai(2*i,k)
1306                ENDDO
1307                DO  i = 1, (nx+1)/2 - 1
1308                   ar(nx+1-i,k) = ai(2*i+1,k)
1309                ENDDO
1310             ENDDO
1311
1312          ELSE
1313
1314             DO  k = 1, nz
1315                DO  i = 0, (nx+1)/2
1316                   ai(2*i,k) = ar(i,k)
1317                ENDDO
1318                DO  i = 1, (nx+1)/2 - 1
1319                   ai(2*i+1,k) = ar(nx+1-i,k)
1320                ENDDO
1321                ai(1,k) = 0.0
1322                ai(nx+2,k) = 0.0
1323             ENDDO
1324
1325             CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
1326
1327             ar(0:nx,1:nz) = ai(0:nx,1:nz)
1328
1329          ENDIF
1330
1331       ELSEIF ( fft_method == 'system-specific' )  THEN
1332
1333#if defined( __nec )
1334          siza = SIZE( ai, 1 )
1335          sizw = SIZE( work, 1 )
1336
1337          IF ( direction == 'forward')  THEN
1338
1339!
1340!--          Tables are initialized once more. This call should not be
1341!--          necessary, but otherwise program aborts in asymmetric case
1342             CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
1343                          trig_xf, work1, 0 )
1344
1345             ai(0:nx,1:nz) = ar(0:nx,1:nz)
1346             IF ( nz1 > nz )  THEN
1347                ai(:,nz1) = 0.0
1348             ENDIF
1349
1350             CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw, &
1351                          trig_xf, work1, 0 )
1352
1353             DO  k = 1, nz
1354                DO  i = 0, (nx+1)/2
1355                   ar(i,k) = REAL( work(i+1,k) )
1356                ENDDO
1357                DO  i = 1, (nx+1)/2 - 1
1358                   ar(nx+1-i,k) = AIMAG( work(i+1,k) )
1359                ENDDO
1360             ENDDO
1361
1362          ELSE
1363
1364!
1365!--          Tables are initialized once more. This call should not be
1366!--          necessary, but otherwise program aborts in asymmetric case
1367             CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
1368                          trig_xb, work1, 0 )
1369
1370             IF ( nz1 > nz )  THEN
1371                work(:,nz1) = 0.0
1372             ENDIF
1373             DO  k = 1, nz
1374                work(1,k) = CMPLX( ar(0,k), 0.0 )
1375                DO  i = 1, (nx+1)/2 - 1
1376                   work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
1377                ENDDO
1378                work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
1379             ENDDO
1380
1381             CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, &
1382                          trig_xb, work1, 0 )
1383
1384             ar(0:nx,1:nz) = ai(0:nx,1:nz)
1385
1386          ENDIF
1387
1388#else
1389          message_string = 'no system-specific fft-call available'
1390          CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 ) 
1391#endif
1392
1393       ELSE
1394
1395          message_string = 'fft method "' // TRIM( fft_method) // &
1396                           '" not available'
1397          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
1398
1399       ENDIF
1400
1401    END SUBROUTINE fft_x_m
1402
1403    SUBROUTINE fft_y_m( ar, ny1, direction )
1404
1405!----------------------------------------------------------------------!
1406!                               fft_y_m                                !
1407!                                                                      !
1408!               Fourier-transformation along y-direction               !
1409!                 Version for 1d domain decomposition                  !
1410!             using multiple 1D FFT from Math Keisan on NEC            !
1411!                       or Temperton-algorithm                         !
1412!       (no singleton-algorithm on NEC because it does not vectorize)  !
1413!                                                                      !
1414!----------------------------------------------------------------------!
1415
1416       IMPLICIT NONE
1417
1418       CHARACTER (LEN=*)              ::  direction
1419       INTEGER                        ::  j, k, ny1, siza
1420
1421       REAL, DIMENSION(0:ny1,nz)      ::  ar
1422       REAL, DIMENSION(0:ny+3,nz+1)   ::  ai
1423       REAL, DIMENSION(6*(ny+4),nz+1) ::  work1
1424#if defined( __nec )
1425       INTEGER                             ::  sizw
1426       COMPLEX, DIMENSION((ny+4)/2+1,nz+1) ::  work
1427#endif
1428
1429       IF ( fft_method == 'temperton-algorithm' )  THEN
1430
1431          siza = SIZE( ai, 1 )
1432
1433          IF ( direction == 'forward')  THEN
1434
1435             ai(0:ny,1:nz) = ar(0:ny,1:nz)
1436             ai(ny+1:,:)   = 0.0
1437
1438             CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
1439
1440             DO  k = 1, nz
1441                DO  j = 0, (ny+1)/2
1442                   ar(j,k) = ai(2*j,k)
1443                ENDDO
1444                DO  j = 1, (ny+1)/2 - 1
1445                   ar(ny+1-j,k) = ai(2*j+1,k)
1446                ENDDO
1447             ENDDO
1448
1449          ELSE
1450
1451             DO  k = 1, nz
1452                DO  j = 0, (ny+1)/2
1453                   ai(2*j,k) = ar(j,k)
1454                ENDDO
1455                DO  j = 1, (ny+1)/2 - 1
1456                   ai(2*j+1,k) = ar(ny+1-j,k)
1457                ENDDO
1458                ai(1,k) = 0.0
1459                ai(ny+2,k) = 0.0
1460             ENDDO
1461
1462             CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
1463
1464             ar(0:ny,1:nz) = ai(0:ny,1:nz)
1465
1466          ENDIF
1467
1468       ELSEIF ( fft_method == 'system-specific' )  THEN
1469
1470#if defined( __nec )
1471          siza = SIZE( ai, 1 )
1472          sizw = SIZE( work, 1 )
1473
1474          IF ( direction == 'forward')  THEN
1475
1476!
1477!--          Tables are initialized once more. This call should not be
1478!--          necessary, but otherwise program aborts in asymmetric case
1479             CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
1480                          trig_yf, work1, 0 )
1481
1482             ai(0:ny,1:nz) = ar(0:ny,1:nz)
1483             IF ( nz1 > nz )  THEN
1484                ai(:,nz1) = 0.0
1485             ENDIF
1486
1487             CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, &
1488                          trig_yf, work1, 0 )
1489
1490             DO  k = 1, nz
1491                DO  j = 0, (ny+1)/2
1492                   ar(j,k) = REAL( work(j+1,k) )
1493                ENDDO
1494                DO  j = 1, (ny+1)/2 - 1
1495                   ar(ny+1-j,k) = AIMAG( work(j+1,k) )
1496                ENDDO
1497             ENDDO
1498
1499          ELSE
1500
1501!
1502!--          Tables are initialized once more. This call should not be
1503!--          necessary, but otherwise program aborts in asymmetric case
1504             CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
1505                          trig_yb, work1, 0 )
1506
1507             IF ( nz1 > nz )  THEN
1508                work(:,nz1) = 0.0
1509             ENDIF
1510             DO  k = 1, nz
1511                work(1,k) = CMPLX( ar(0,k), 0.0 )
1512                DO  j = 1, (ny+1)/2 - 1
1513                   work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
1514                ENDDO
1515                work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
1516             ENDDO
1517
1518             CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, &
1519                          trig_yb, work1, 0 )
1520
1521             ar(0:ny,1:nz) = ai(0:ny,1:nz)
1522
1523          ENDIF
1524
1525#else
1526          message_string = 'no system-specific fft-call available'
1527          CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 ) 
1528#endif
1529
1530       ELSE
1531         
1532          message_string = 'fft method "' // TRIM( fft_method) // &
1533                           '" not available'
1534          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
1535
1536       ENDIF
1537
1538    END SUBROUTINE fft_y_m
1539
1540
1541 END MODULE fft_xy
Note: See TracBrowser for help on using the repository browser.