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

Last change on this file since 1115 was 1112, checked in by raasch, 12 years ago

last commit documented

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