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

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

last commit documented

  • Property svn:keywords set to Id
File size: 25.5 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 1037 2012-10-22 14:10:22Z raasch $
27!
28! 1036 2012-10-22 13:43:42Z raasch
29! code put under GPL (PALM 3.9)
30!
31! 274 2009-03-26 15:11:21Z heinze
32! Output of messages replaced by message handling routine.
33!
34! Feb. 2007
35! RCS Log replace by Id keyword, revision history cleaned up
36!
37! Revision 1.4  2006/03/28 12:27:09  raasch
38! Stop when system-specific fft is selected on NEC. For unknown reasons this
39! causes a program abort during first allocation in init_grid.
40!
41! Revision 1.2  2004/04/30 11:44:27  raasch
42! Module renamed from fft_for_1d_decomp to fft_xy, 1d-routines renamed to
43! fft_x and fft_y,
44! function FFT replaced by subroutine FFTN due to problems with 64-bit
45! mode on ibm,
46! shape of array cwork is explicitly stored in ishape/jshape and handled
47! to routine FFTN instead of shape-function (due to compiler error on
48! decalpha),
49! non vectorized FFT for nec included
50!
51! Revision 1.1  2002/06/11 13:00:49  raasch
52! Initial revision
53!
54!
55! Description:
56! ------------
57! Fast Fourier transformation along x and y for 1d domain decomposition along x.
58! Original version: Klaus Ketelsen (May 2002)
59!------------------------------------------------------------------------------!
60
61    USE array_kind
62    USE control_parameters
63    USE indices
64    USE singleton
65    USE temperton_fft
66
67    IMPLICIT NONE
68
69    PRIVATE
70    PUBLIC fft_x, fft_y, fft_init, fft_x_m, fft_y_m
71
72    INTEGER, DIMENSION(:), ALLOCATABLE, SAVE ::  ifax_x, ifax_y
73
74    LOGICAL, SAVE                            ::  init_fft = .FALSE.
75
76    REAL, SAVE ::  sqr_nx, sqr_ny
77    REAL, DIMENSION(:), ALLOCATABLE, SAVE    ::  trigs_x, trigs_y
78
79#if defined( __ibm )
80    INTEGER, PARAMETER ::  nau1 = 20000, nau2 = 22000
81!
82!-- The following working arrays contain tables and have to be "save" and
83!-- shared in OpenMP sense
84    REAL, DIMENSION(nau1), SAVE ::  aux1, auy1, aux3, auy3
85#elif defined( __nec )
86    INTEGER, SAVE ::  nz1
87    REAL, DIMENSION(:), ALLOCATABLE, SAVE ::  trig_xb, trig_xf, trig_yb, &
88                                              trig_yf
89#endif
90
91!
92!-- Public interfaces
93    INTERFACE fft_init
94       MODULE PROCEDURE fft_init
95    END INTERFACE fft_init
96
97    INTERFACE fft_x
98       MODULE PROCEDURE fft_x
99    END INTERFACE fft_x
100
101    INTERFACE fft_y
102       MODULE PROCEDURE fft_y
103    END INTERFACE fft_y
104
105    INTERFACE fft_x_m
106       MODULE PROCEDURE fft_x_m
107    END INTERFACE fft_x_m
108
109    INTERFACE fft_y_m
110       MODULE PROCEDURE fft_y_m
111    END INTERFACE fft_y_m
112
113 CONTAINS
114
115
116    SUBROUTINE fft_init
117
118       IMPLICIT NONE
119
120!
121!--    The following temporary working arrays have to be on stack or private
122!--    in OpenMP sense
123#if defined( __ibm )
124       REAL, DIMENSION(0:nx+2) :: workx
125       REAL, DIMENSION(0:ny+2) :: worky
126       REAL, DIMENSION(nau2)   :: aux2, auy2, aux4, auy4
127#elif defined( __nec )
128       REAL, DIMENSION(0:nx+3,nz+1)   ::  work_x
129       REAL, DIMENSION(0:ny+3,nz+1)   ::  work_y
130       REAL, DIMENSION(6*(nx+3),nz+1) ::  workx
131       REAL, DIMENSION(6*(ny+3),nz+1) ::  worky
132#endif
133
134!
135!--    Return, if already called
136       IF ( init_fft )  THEN
137          RETURN
138       ELSE
139          init_fft = .TRUE.
140       ENDIF
141
142       IF ( fft_method == 'system-specific' )  THEN
143
144          sqr_nx = SQRT( 1.0 / ( nx + 1.0 ) )
145          sqr_ny = SQRT( 1.0 / ( ny + 1.0 ) )
146#if defined( __ibm )  &&  ! defined( __ibmy_special )
147!
148!--       Initialize tables for fft along x
149          CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1,  1, sqr_nx, aux1, nau1, &
150                      aux2, nau2 )
151          CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_nx, aux3, nau1, &
152                      aux4, nau2 )
153!
154!--       Initialize tables for fft along y
155          CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1,  1, sqr_ny, auy1, nau1, &
156                      auy2, nau2 )
157          CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_ny, auy3, nau1, &
158                      auy4, nau2 )
159#elif defined( __nec )
160          message_string = 'fft method "' // TRIM( fft_method) // &
161                           '" currently does not work on NEC'
162          CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 )
163
164          ALLOCATE( trig_xb(2*(nx+1)), trig_xf(2*(nx+1)), &
165                    trig_yb(2*(ny+1)), trig_yf(2*(ny+1)) )
166
167          work_x = 0.0
168          work_y = 0.0
169          nz1  = nz + MOD( nz+1, 2 )  ! odd nz slows down fft significantly
170                                      ! when using the NEC ffts
171
172!
173!--       Initialize tables for fft along x (non-vector and vector case (M))
174          CALL DZFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xf, workx, 0 )
175          CALL ZDFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xb, workx, 0 )
176          CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
177                       trig_xf, workx, 0 )
178          CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
179                       trig_xb, workx, 0 )
180!
181!--       Initialize tables for fft along y (non-vector and vector case (M))
182          CALL DZFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yf, worky, 0 )
183          CALL ZDFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yb, worky, 0 )
184          CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
185                       trig_yf, worky, 0 )
186          CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
187                       trig_yb, worky, 0 )
188#else
189          message_string = 'no system-specific fft-call available'
190          CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
191#endif
192       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
193!
194!--       Temperton-algorithm
195!--       Initialize tables for fft along x and y
196          ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
197
198          CALL set99( trigs_x, ifax_x, nx+1 )
199          CALL set99( trigs_y, ifax_y, ny+1 )
200
201       ELSEIF ( fft_method == 'singleton-algorithm' )  THEN
202
203          CONTINUE
204
205       ELSE
206
207          message_string = 'fft method "' // TRIM( fft_method) // &
208                           '" not available'
209          CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
210       ENDIF
211
212    END SUBROUTINE fft_init
213
214
215    SUBROUTINE fft_x( ar, direction )
216
217!----------------------------------------------------------------------!
218!                                 fft_x                                !
219!                                                                      !
220!               Fourier-transformation along x-direction               !
221!                                                                      !
222!      fft_x uses internal algorithms (Singleton or Temperton) or      !
223!           system-specific routines, if they are available            !
224!----------------------------------------------------------------------!
225
226       IMPLICIT NONE
227
228       CHARACTER (LEN=*) ::  direction
229       INTEGER ::  i, ishape(1)
230
231!kk    REAL, DIMENSION(:)        ::  ar !kk Does NOT work (Bug??)
232       REAL, DIMENSION(0:nx)     ::  ar
233       REAL, DIMENSION(0:nx+2)   ::  work
234       REAL, DIMENSION(nx+2)     ::  work1
235       COMPLEX, DIMENSION(:), ALLOCATABLE ::  cwork
236#if defined( __ibm )
237       REAL, DIMENSION(nau2)     ::  aux2, aux4
238#elif defined( __nec )
239       REAL, DIMENSION(6*(nx+1)) ::  work2
240#endif
241
242       IF ( fft_method == 'singleton-algorithm' )  THEN
243
244!
245!--       Performing the fft with singleton's software works on every system,
246!--       since it is part of the model
247          ALLOCATE( cwork(0:nx) )
248     
249          IF ( direction == 'forward')   then
250
251             DO  i = 0, nx
252                cwork(i) = CMPLX( ar(i) )
253             ENDDO
254             ishape = SHAPE( cwork )
255             CALL FFTN( cwork, ishape )
256
257             DO  i = 0, (nx+1)/2
258                ar(i) = REAL( cwork(i) )
259             ENDDO
260             DO  i = 1, (nx+1)/2 - 1
261                ar(nx+1-i) = -AIMAG( cwork(i) )
262             ENDDO
263
264          ELSE
265
266             cwork(0) = CMPLX( ar(0), 0.0 )
267             DO  i = 1, (nx+1)/2 - 1
268                cwork(i)      = CMPLX( ar(i), -ar(nx+1-i) )
269                cwork(nx+1-i) = CMPLX( ar(i),  ar(nx+1-i) )
270             ENDDO
271             cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
272
273             ishape = SHAPE( cwork )
274             CALL FFTN( cwork, ishape, inv = .TRUE. )
275
276             DO  i = 0, nx
277                ar(i) = REAL( cwork(i) )
278             ENDDO
279
280          ENDIF
281
282          DEALLOCATE( cwork )
283
284       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
285
286!
287!--       Performing the fft with Temperton's software works on every system,
288!--       since it is part of the model
289          IF ( direction == 'forward' )  THEN
290
291             work(0:nx) = ar
292             CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
293
294             DO  i = 0, (nx+1)/2
295                ar(i) = work(2*i)
296             ENDDO
297             DO  i = 1, (nx+1)/2 - 1
298                ar(nx+1-i) = work(2*i+1)
299             ENDDO
300
301          ELSE
302
303             DO  i = 0, (nx+1)/2
304                work(2*i) = ar(i)
305             ENDDO
306             DO  i = 1, (nx+1)/2 - 1
307                work(2*i+1) = ar(nx+1-i)
308             ENDDO
309             work(1)    = 0.0
310             work(nx+2) = 0.0
311
312             CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
313             ar = work(0:nx)
314
315          ENDIF
316
317       ELSEIF ( fft_method == 'system-specific' )  THEN
318
319#if defined( __ibm )  &&  ! defined( __ibmy_special )
320          IF ( direction == 'forward' )  THEN
321
322             CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_nx, aux1, nau1, &
323                         aux2, nau2 )
324
325             DO  i = 0, (nx+1)/2
326                ar(i) = work(2*i)
327             ENDDO
328             DO  i = 1, (nx+1)/2 - 1
329                ar(nx+1-i) = work(2*i+1)
330             ENDDO
331
332          ELSE
333
334             DO  i = 0, (nx+1)/2
335                work(2*i) = ar(i)
336             ENDDO
337             DO  i = 1, (nx+1)/2 - 1
338                work(2*i+1) = ar(nx+1-i)
339             ENDDO
340             work(1) = 0.0
341             work(nx+2) = 0.0
342
343             CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_nx, aux3, nau1, &
344                         aux4, nau2 )
345
346             DO  i = 0, nx
347                ar(i) = work(i)
348             ENDDO
349
350          ENDIF
351#elif defined( __nec )
352          IF ( direction == 'forward' )  THEN
353
354             work(0:nx) = ar(0:nx)
355
356             CALL DZFFT( 1, nx+1, sqr_nx, work, work, trig_xf, work2, 0 )
357
358             DO  i = 0, (nx+1)/2
359                ar(i) = work(2*i)
360             ENDDO
361             DO  i = 1, (nx+1)/2 - 1
362                ar(nx+1-i) = work(2*i+1)
363             ENDDO
364
365          ELSE
366
367             DO  i = 0, (nx+1)/2
368                work(2*i) = ar(i)
369             ENDDO
370             DO  i = 1, (nx+1)/2 - 1
371                work(2*i+1) = ar(nx+1-i)
372             ENDDO
373             work(1) = 0.0
374             work(nx+2) = 0.0
375
376             CALL ZDFFT( -1, nx+1, sqr_nx, work, work, trig_xb, work2, 0 )
377
378             ar(0:nx) = work(0:nx)
379
380          ENDIF
381#else
382          message_string = 'no system-specific fft-call available'
383          CALL message( 'fft_x', 'PA0188', 1, 2, 0, 6, 0 )
384#endif
385       ELSE
386          message_string = 'fft method "' // TRIM( fft_method) // &
387                           '" not available'
388          CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
389
390       ENDIF
391
392    END SUBROUTINE fft_x
393
394    SUBROUTINE fft_y( ar, direction )
395
396!----------------------------------------------------------------------!
397!                                 fft_y                                !
398!                                                                      !
399!               Fourier-transformation along y-direction               !
400!                                                                      !
401!      fft_y uses internal algorithms (Singleton or Temperton) or      !
402!           system-specific routines, if they are available            !
403!----------------------------------------------------------------------!
404
405       IMPLICIT NONE
406
407       CHARACTER (LEN=*) ::  direction
408       INTEGER ::  j, jshape(1)
409
410!kk    REAL, DIMENSION(:)        ::  ar !kk Does NOT work (Bug??)
411       REAL, DIMENSION(0:ny)     ::  ar
412       REAL, DIMENSION(0:ny+2)   ::  work
413       REAL, DIMENSION(ny+2)     ::  work1
414       COMPLEX, DIMENSION(:), ALLOCATABLE ::  cwork
415#if defined( __ibm )
416       REAL, DIMENSION(nau2)     ::  auy2, auy4
417#elif defined( __nec )
418       REAL, DIMENSION(6*(ny+1)) ::  work2
419#endif
420
421       IF ( fft_method == 'singleton-algorithm' )  THEN
422
423!
424!--       Performing the fft with singleton's software works on every system,
425!--       since it is part of the model
426          ALLOCATE( cwork(0:ny) )
427
428          IF ( direction == 'forward')  THEN
429
430             DO  j = 0, ny
431                cwork(j) = CMPLX( ar(j) )
432             ENDDO
433
434             jshape = SHAPE( cwork )
435             CALL FFTN( cwork, jshape )
436
437             DO  j = 0, (ny+1)/2
438                ar(j) = REAL( cwork(j) )
439             ENDDO
440             DO  j = 1, (ny+1)/2 - 1
441                ar(ny+1-j) = -AIMAG( cwork(j) )
442             ENDDO
443
444          ELSE
445
446             cwork(0) = CMPLX( ar(0), 0.0 )
447             DO  j = 1, (ny+1)/2 - 1
448                cwork(j)      = CMPLX( ar(j), -ar(ny+1-j) )
449                cwork(ny+1-j) = CMPLX( ar(j),  ar(ny+1-j) )
450             ENDDO
451             cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
452
453             jshape = SHAPE( cwork )
454             CALL FFTN( cwork, jshape, inv = .TRUE. )
455
456             DO  j = 0, ny
457                ar(j) = REAL( cwork(j) )
458             ENDDO
459
460          ENDIF
461
462          DEALLOCATE( cwork )
463
464       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
465
466!
467!--       Performing the fft with Temperton's software works on every system,
468!--       since it is part of the model
469          IF ( direction == 'forward' )  THEN
470
471             work(0:ny) = ar
472             CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
473
474             DO  j = 0, (ny+1)/2
475                ar(j) = work(2*j)
476             ENDDO
477             DO  j = 1, (ny+1)/2 - 1
478                ar(ny+1-j) = work(2*j+1)
479             ENDDO
480
481          ELSE
482
483             DO  j = 0, (ny+1)/2
484                work(2*j) = ar(j)
485             ENDDO
486             DO  j = 1, (ny+1)/2 - 1
487                work(2*j+1) = ar(ny+1-j)
488             ENDDO
489             work(1)    = 0.0
490             work(ny+2) = 0.0
491
492             CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
493             ar = work(0:ny)
494
495          ENDIF
496
497       ELSEIF ( fft_method == 'system-specific' )  THEN
498
499#if defined( __ibm )  &&  ! defined( __ibmy_special )
500          IF ( direction == 'forward')  THEN
501
502             CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_ny, auy1, nau1, &
503                         auy2, nau2 )
504
505             DO  j = 0, (ny+1)/2
506                ar(j) = work(2*j)
507             ENDDO
508             DO  j = 1, (ny+1)/2 - 1
509                ar(ny+1-j) = work(2*j+1)
510             ENDDO
511
512          ELSE
513
514             DO  j = 0, (ny+1)/2
515                work(2*j) = ar(j)
516             ENDDO
517             DO  j = 1, (ny+1)/2 - 1
518                work(2*j+1) = ar(ny+1-j)
519             ENDDO
520             work(1)    = 0.0
521             work(ny+2) = 0.0
522
523             CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_ny, auy3, nau1, &
524                         auy4, nau2 )
525
526             DO  j = 0, ny
527                ar(j) = work(j)
528             ENDDO
529
530          ENDIF
531#elif defined( __nec )
532          IF ( direction == 'forward' )  THEN
533
534             work(0:ny) = ar(0:ny)
535
536             CALL DZFFT( 1, ny+1, sqr_ny, work, work, trig_yf, work2, 0 )
537
538             DO  j = 0, (ny+1)/2
539                ar(j) = work(2*j)
540             ENDDO
541             DO  j = 1, (ny+1)/2 - 1
542                ar(ny+1-j) = work(2*j+1)
543             ENDDO
544
545          ELSE
546
547             DO  j = 0, (ny+1)/2
548                work(2*j) = ar(j)
549             ENDDO
550             DO  j = 1, (ny+1)/2 - 1
551                work(2*j+1) = ar(ny+1-j)
552             ENDDO
553             work(1) = 0.0
554             work(ny+2) = 0.0
555
556             CALL ZDFFT( -1, ny+1, sqr_ny, work, work, trig_yb, work2, 0 )
557
558             ar(0:ny) = work(0:ny)
559
560          ENDIF
561#else
562          message_string = 'no system-specific fft-call available'
563          CALL message( 'fft_y', 'PA0188', 1, 2, 0, 6, 0 ) 
564
565#endif
566
567       ELSE
568
569          message_string = 'fft method "' // TRIM( fft_method) // &
570                           '" not available'
571          CALL message( 'fft_y', 'PA0189', 1, 2, 0, 6, 0 )
572
573       ENDIF
574
575    END SUBROUTINE fft_y
576
577    SUBROUTINE fft_x_m( ar, direction )
578
579!----------------------------------------------------------------------!
580!                               fft_x_m                                !
581!                                                                      !
582!               Fourier-transformation along x-direction               !
583!                 Version for 1d domain decomposition                  !
584!             using multiple 1D FFT from Math Keisan on NEC            !
585!                       or Temperton-algorithm                         !
586!       (no singleton-algorithm on NEC because it does not vectorize)  !
587!                                                                      !
588!----------------------------------------------------------------------!
589
590       IMPLICIT NONE
591
592       CHARACTER (LEN=*)              ::  direction
593       INTEGER                        ::  i, k, siza, sizw
594
595       REAL, DIMENSION(0:nx,nz)       ::  ar
596       REAL, DIMENSION(0:nx+3,nz+1)   ::  ai
597       REAL, DIMENSION(6*(nx+4),nz+1) ::  work1
598#if defined( __nec )
599       COMPLEX, DIMENSION((nx+4)/2+1,nz+1) ::  work
600#endif
601
602       IF ( fft_method == 'temperton-algorithm' )  THEN
603
604          siza = SIZE( ai, 1 )
605
606          IF ( direction == 'forward')  THEN
607
608             ai(0:nx,1:nz) = ar(0:nx,1:nz)
609             ai(nx+1:,:)   = 0.0
610
611             CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
612
613             DO  k = 1, nz
614                DO  i = 0, (nx+1)/2
615                   ar(i,k) = ai(2*i,k)
616                ENDDO
617                DO  i = 1, (nx+1)/2 - 1
618                   ar(nx+1-i,k) = ai(2*i+1,k)
619                ENDDO
620             ENDDO
621
622          ELSE
623
624             DO  k = 1, nz
625                DO  i = 0, (nx+1)/2
626                   ai(2*i,k) = ar(i,k)
627                ENDDO
628                DO  i = 1, (nx+1)/2 - 1
629                   ai(2*i+1,k) = ar(nx+1-i,k)
630                ENDDO
631                ai(1,k) = 0.0
632                ai(nx+2,k) = 0.0
633             ENDDO
634
635             CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
636
637             ar(0:nx,1:nz) = ai(0:nx,1:nz)
638
639          ENDIF
640
641       ELSEIF ( fft_method == 'system-specific' )  THEN
642
643#if defined( __nec )
644          siza = SIZE( ai, 1 )
645          sizw = SIZE( work, 1 )
646
647          IF ( direction == 'forward')  THEN
648
649!
650!--          Tables are initialized once more. This call should not be
651!--          necessary, but otherwise program aborts in asymmetric case
652             CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
653                          trig_xf, work1, 0 )
654
655             ai(0:nx,1:nz) = ar(0:nx,1:nz)
656             IF ( nz1 > nz )  THEN
657                ai(:,nz1) = 0.0
658             ENDIF
659
660             CALL DZFFTM( 1, nx+1, nz1, sqr_nx, ai, siza, work, sizw, &
661                          trig_xf, work1, 0 )
662
663             DO  k = 1, nz
664                DO  i = 0, (nx+1)/2
665                   ar(i,k) = REAL( work(i+1,k) )
666                ENDDO
667                DO  i = 1, (nx+1)/2 - 1
668                   ar(nx+1-i,k) = AIMAG( work(i+1,k) )
669                ENDDO
670             ENDDO
671
672          ELSE
673
674!
675!--          Tables are initialized once more. This call should not be
676!--          necessary, but otherwise program aborts in asymmetric case
677             CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
678                          trig_xb, work1, 0 )
679
680             IF ( nz1 > nz )  THEN
681                work(:,nz1) = 0.0
682             ENDIF
683             DO  k = 1, nz
684                work(1,k) = CMPLX( ar(0,k), 0.0 )
685                DO  i = 1, (nx+1)/2 - 1
686                   work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
687                ENDDO
688                work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
689             ENDDO
690
691             CALL ZDFFTM( -1, nx+1, nz1, sqr_nx, work, sizw, ai, siza, &
692                          trig_xb, work1, 0 )
693
694             ar(0:nx,1:nz) = ai(0:nx,1:nz)
695
696          ENDIF
697
698#else
699          message_string = 'no system-specific fft-call available'
700          CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 ) 
701#endif
702
703       ELSE
704
705          message_string = 'fft method "' // TRIM( fft_method) // &
706                           '" not available'
707          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
708
709       ENDIF
710
711    END SUBROUTINE fft_x_m
712
713    SUBROUTINE fft_y_m( ar, ny1, direction )
714
715!----------------------------------------------------------------------!
716!                               fft_y_m                                !
717!                                                                      !
718!               Fourier-transformation along y-direction               !
719!                 Version for 1d domain decomposition                  !
720!             using multiple 1D FFT from Math Keisan on NEC            !
721!                       or Temperton-algorithm                         !
722!       (no singleton-algorithm on NEC because it does not vectorize)  !
723!                                                                      !
724!----------------------------------------------------------------------!
725
726       IMPLICIT NONE
727
728       CHARACTER (LEN=*)              ::  direction
729       INTEGER                        ::  j, k, ny1, siza, sizw
730
731       REAL, DIMENSION(0:ny1,nz)      ::  ar
732       REAL, DIMENSION(0:ny+3,nz+1)   ::  ai
733       REAL, DIMENSION(6*(ny+4),nz+1) ::  work1
734#if defined( __nec )
735       COMPLEX, DIMENSION((ny+4)/2+1,nz+1) ::  work
736#endif
737
738       IF ( fft_method == 'temperton-algorithm' )  THEN
739
740          siza = SIZE( ai, 1 )
741
742          IF ( direction == 'forward')  THEN
743
744             ai(0:ny,1:nz) = ar(0:ny,1:nz)
745             ai(ny+1:,:)   = 0.0
746
747             CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
748
749             DO  k = 1, nz
750                DO  j = 0, (ny+1)/2
751                   ar(j,k) = ai(2*j,k)
752                ENDDO
753                DO  j = 1, (ny+1)/2 - 1
754                   ar(ny+1-j,k) = ai(2*j+1,k)
755                ENDDO
756             ENDDO
757
758          ELSE
759
760             DO  k = 1, nz
761                DO  j = 0, (ny+1)/2
762                   ai(2*j,k) = ar(j,k)
763                ENDDO
764                DO  j = 1, (ny+1)/2 - 1
765                   ai(2*j+1,k) = ar(ny+1-j,k)
766                ENDDO
767                ai(1,k) = 0.0
768                ai(ny+2,k) = 0.0
769             ENDDO
770
771             CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
772
773             ar(0:ny,1:nz) = ai(0:ny,1:nz)
774
775          ENDIF
776
777       ELSEIF ( fft_method == 'system-specific' )  THEN
778
779#if defined( __nec )
780          siza = SIZE( ai, 1 )
781          sizw = SIZE( work, 1 )
782
783          IF ( direction == 'forward')  THEN
784
785!
786!--          Tables are initialized once more. This call should not be
787!--          necessary, but otherwise program aborts in asymmetric case
788             CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
789                          trig_yf, work1, 0 )
790
791             ai(0:ny,1:nz) = ar(0:ny,1:nz)
792             IF ( nz1 > nz )  THEN
793                ai(:,nz1) = 0.0
794             ENDIF
795
796             CALL DZFFTM( 1, ny+1, nz1, sqr_ny, ai, siza, work, sizw, &
797                          trig_yf, work1, 0 )
798
799             DO  k = 1, nz
800                DO  j = 0, (ny+1)/2
801                   ar(j,k) = REAL( work(j+1,k) )
802                ENDDO
803                DO  j = 1, (ny+1)/2 - 1
804                   ar(ny+1-j,k) = AIMAG( work(j+1,k) )
805                ENDDO
806             ENDDO
807
808          ELSE
809
810!
811!--          Tables are initialized once more. This call should not be
812!--          necessary, but otherwise program aborts in asymmetric case
813             CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
814                          trig_yb, work1, 0 )
815
816             IF ( nz1 > nz )  THEN
817                work(:,nz1) = 0.0
818             ENDIF
819             DO  k = 1, nz
820                work(1,k) = CMPLX( ar(0,k), 0.0 )
821                DO  j = 1, (ny+1)/2 - 1
822                   work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
823                ENDDO
824                work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
825             ENDDO
826
827             CALL ZDFFTM( -1, ny+1, nz1, sqr_ny, work, sizw, ai, siza, &
828                          trig_yb, work1, 0 )
829
830             ar(0:ny,1:nz) = ai(0:ny,1:nz)
831
832          ENDIF
833
834#else
835          message_string = 'no system-specific fft-call available'
836          CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 ) 
837#endif
838
839       ELSE
840         
841          message_string = 'fft method "' // TRIM( fft_method) // &
842                           '" not available'
843          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
844
845       ENDIF
846
847    END SUBROUTINE fft_y_m
848
849 END MODULE fft_xy
Note: See TracBrowser for help on using the repository browser.