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

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

last commit documented

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