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

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