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

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

last commit documented

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