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

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

code adjustment of accelerator version for PGI 12.3 / CUDA 5.0

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