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

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

last commit documented

  • Property svn:keywords set to Id
File size: 45.7 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!
23!
24! Former revisions:
25! -----------------
26! $Id: fft_xy.f90 1112 2013-03-09 00:34:37Z raasch $
27!
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!
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!
38! 1092 2013-02-02 11:24:22Z raasch
39! variable sizw declared for NEC case only
40!
41! 1036 2012-10-22 13:43:42Z raasch
42! code put under GPL (PALM 3.9)
43!
44! 274 2009-03-26 15:11:21Z heinze
45! Output of messages replaced by message handling routine.
46!
47! Feb. 2007
48! RCS Log replace by Id keyword, revision history cleaned up
49!
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
76    USE precision_kind
77    USE singleton
78    USE temperton_fft
79    USE transpose_indices
80
81    IMPLICIT NONE
82
83    PRIVATE
84    PUBLIC fft_x, fft_x_1d, fft_y, fft_y_1d, fft_init, fft_x_m, fft_y_m
85
86    INTEGER, DIMENSION(:), ALLOCATABLE, SAVE ::  ifax_x, ifax_y
87
88    LOGICAL, SAVE                            ::  init_fft = .FALSE.
89
90    REAL, SAVE ::  dnx, dny, sqr_dnx, sqr_dny
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
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
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
118    INTERFACE fft_x_1d
119       MODULE PROCEDURE fft_x_1d
120    END INTERFACE fft_x_1d
121
122    INTERFACE fft_y
123       MODULE PROCEDURE fft_y
124    END INTERFACE fft_y
125
126    INTERFACE fft_y_1d
127       MODULE PROCEDURE fft_y_1d
128    END INTERFACE fft_y_1d
129
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
143       USE cuda_fft_interfaces
144
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
171          dnx = 1.0 / ( nx + 1.0 )
172          dny = 1.0 / ( ny + 1.0 )
173          sqr_dnx = SQRT( dnx )
174          sqr_dny = SQRT( dny )
175#if defined( __ibm )  &&  ! defined( __ibmy_special )
176!
177!--       Initialize tables for fft along x
178          CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1,  1, sqr_dnx, aux1, nau1, &
179                      aux2, nau2 )
180          CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
181                      aux4, nau2 )
182!
183!--       Initialize tables for fft along y
184          CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1,  1, sqr_dny, auy1, nau1, &
185                      auy2, nau2 )
186          CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
187                      auy4, nau2 )
188#elif defined( __nec )
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 )
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))
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, &
206                       trig_xf, workx, 0 )
207          CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
208                       trig_xb, workx, 0 )
209!
210!--       Initialize tables for fft along y (non-vector and vector case (M))
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, &
214                       trig_yf, worky, 0 )
215          CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
216                       trig_yb, worky, 0 )
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)
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) )
224#else
225          message_string = 'no system-specific fft-call available'
226          CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
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
243          message_string = 'fft method "' // TRIM( fft_method) // &
244                           '" not available'
245          CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
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               !
257!                     Version for 2D-decomposition                     !
258!                                                                      !
259!      fft_x uses internal algorithms (Singleton or Temperton) or      !
260!           system-specific routines, if they are available            !
261!----------------------------------------------------------------------!
262
263       USE cuda_fft_interfaces
264
265       IMPLICIT NONE
266
267       CHARACTER (LEN=*) ::  direction
268       INTEGER ::  i, ishape(1), j, k
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 )
280       !$acc declare create( ar_tmp )
281       COMPLEX(dpk), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) ::  ar_tmp
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
510             !$acc data present( ar )
511             CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
512
513             !$acc kernels
514             !$acc loop
515             DO  k = nzb_x, nzt_x
516                DO  j = nys_x, nyn_x
517
518                   !$acc loop vector( 32 )
519                   DO  i = 0, (nx+1)/2
520                      ar(i,j,k)      = REAL( ar_tmp(i,j,k) )  * dnx
521                   ENDDO
522
523                   !$acc loop vector( 32 )
524                   DO  i = 1, (nx+1)/2 - 1
525                      ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
526                   ENDDO
527
528                ENDDO
529             ENDDO
530             !$acc end kernels
531             !$acc end data
532
533          ELSE
534
535             !$acc data present( ar )
536             !$acc kernels
537             !$acc loop
538             DO  k = nzb_x, nzt_x
539                DO  j = nys_x, nyn_x
540
541                   ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
542
543                   !$acc loop vector( 32 )
544                   DO  i = 1, (nx+1)/2 - 1
545                      ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
546                   ENDDO
547                   ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
548
549                ENDDO
550             ENDDO
551             !$acc end kernels
552
553             CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
554             !$acc end data
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
588       INTEGER ::  i, ishape(1)
589
590       LOGICAL ::  forward_fft
591
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
602       IF ( direction == 'forward' )  THEN
603          forward_fft = .TRUE.
604       ELSE
605          forward_fft = .FALSE.
606       ENDIF
607
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     
615          IF ( forward_fft )   then
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
654          IF ( forward_fft )  THEN
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 )
685          IF ( forward_fft )  THEN
686
687             CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
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
708             CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
709                         aux4, nau2 )
710
711             DO  i = 0, nx
712                ar(i) = work(i)
713             ENDDO
714
715          ENDIF
716#elif defined( __nec )
717          IF ( forward_fft )  THEN
718
719             work(0:nx) = ar(0:nx)
720
721             CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
722     
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
741             CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
742
743             ar(0:nx) = work(0:nx)
744
745          ENDIF
746#else
747          message_string = 'no system-specific fft-call available'
748          CALL message( 'fft_x_1d', 'PA0188', 1, 2, 0, 6, 0 )
749#endif
750       ELSE
751          message_string = 'fft method "' // TRIM( fft_method) // &
752                           '" not available'
753          CALL message( 'fft_x_1d', 'PA0189', 1, 2, 0, 6, 0 )
754
755       ENDIF
756
757    END SUBROUTINE fft_x_1d
758
759    SUBROUTINE fft_y( ar, direction )
760
761!----------------------------------------------------------------------!
762!                                 fft_y                                !
763!                                                                      !
764!               Fourier-transformation along y-direction               !
765!                     Version for 2D-decomposition                     !
766!                                                                      !
767!      fft_y uses internal algorithms (Singleton or Temperton) or      !
768!           system-specific routines, if they are available            !
769!----------------------------------------------------------------------!
770
771       USE cuda_fft_interfaces
772
773       IMPLICIT NONE
774
775       CHARACTER (LEN=*) ::  direction
776       INTEGER ::  i, j, jshape(1), k
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 )
788       !$acc declare create( ar_tmp )
789       COMPLEX(dpk), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) ::  ar_tmp
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
1015             !$acc data present( ar )
1016             CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
1017
1018             !$acc kernels
1019             !$acc loop
1020             DO  k = nzb_y, nzt_y
1021                DO  i = nxl_y, nxr_y
1022
1023                   !$acc loop vector( 32 )
1024                   DO  j = 0, (ny+1)/2
1025                      ar(j,i,k)      = REAL( ar_tmp(j,i,k) )  * dny
1026                   ENDDO
1027
1028                   !$acc loop vector( 32 )
1029                   DO  j = 1, (ny+1)/2 - 1
1030                      ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
1031                   ENDDO
1032
1033                ENDDO
1034             ENDDO
1035             !$acc end kernels
1036             !$acc end data
1037
1038          ELSE
1039
1040             !$acc data present( ar )
1041             !$acc kernels
1042             !$acc loop
1043             DO  k = nzb_y, nzt_y
1044                DO  i = nxl_y, nxr_y
1045
1046                   ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
1047
1048                   !$acc loop vector( 32 )
1049                   DO  j = 1, (ny+1)/2 - 1
1050                      ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
1051                   ENDDO
1052                   ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
1053
1054                ENDDO
1055             ENDDO
1056             !$acc end kernels
1057
1058             CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
1059             !$acc end data
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
1093       INTEGER ::  j, jshape(1)
1094
1095       LOGICAL ::  forward_fft
1096
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
1107       IF ( direction == 'forward' )  THEN
1108          forward_fft = .TRUE.
1109       ELSE
1110          forward_fft = .FALSE.
1111       ENDIF
1112
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
1120          IF ( forward_fft )  THEN
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
1161          IF ( forward_fft )  THEN
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 )
1192          IF ( forward_fft )  THEN
1193
1194             CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
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
1215             CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
1216                         auy4, nau2 )
1217
1218             DO  j = 0, ny
1219                ar(j) = work(j)
1220             ENDDO
1221
1222          ENDIF
1223#elif defined( __nec )
1224          IF ( forward_fft )  THEN
1225
1226             work(0:ny) = ar(0:ny)
1227
1228             CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
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
1248             CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1249
1250             ar(0:ny) = work(0:ny)
1251
1252          ENDIF
1253#else
1254          message_string = 'no system-specific fft-call available'
1255          CALL message( 'fft_y_1d', 'PA0188', 1, 2, 0, 6, 0 ) 
1256
1257#endif
1258
1259       ELSE
1260
1261          message_string = 'fft method "' // TRIM( fft_method) // &
1262                           '" not available'
1263          CALL message( 'fft_y_1d', 'PA0189', 1, 2, 0, 6, 0 )
1264
1265       ENDIF
1266
1267    END SUBROUTINE fft_y_1d
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
1285       INTEGER                        ::  i, k, siza
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 )
1291       INTEGER                             ::  sizw
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
1345             CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
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
1353             CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw, &
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
1370             CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
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
1384             CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, &
1385                          trig_xb, work1, 0 )
1386
1387             ar(0:nx,1:nz) = ai(0:nx,1:nz)
1388
1389          ENDIF
1390
1391#else
1392          message_string = 'no system-specific fft-call available'
1393          CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 ) 
1394#endif
1395
1396       ELSE
1397
1398          message_string = 'fft method "' // TRIM( fft_method) // &
1399                           '" not available'
1400          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
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
1422       INTEGER                        ::  j, k, ny1, siza
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 )
1428       INTEGER                             ::  sizw
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
1482             CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
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
1490             CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, &
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
1507             CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
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
1521             CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, &
1522                          trig_yb, work1, 0 )
1523
1524             ar(0:ny,1:nz) = ai(0:ny,1:nz)
1525
1526          ENDIF
1527
1528#else
1529          message_string = 'no system-specific fft-call available'
1530          CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 ) 
1531#endif
1532
1533       ELSE
1534         
1535          message_string = 'fft method "' // TRIM( fft_method) // &
1536                           '" not available'
1537          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
1538
1539       ENDIF
1540
1541    END SUBROUTINE fft_y_m
1542
1543
1544 END MODULE fft_xy
Note: See TracBrowser for help on using the repository browser.