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

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

unused variables remove from several routines

  • Property svn:keywords set to Id
File size: 25.7 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! variable sizw declared for NEC case only
23!
24! Former revisions:
25! -----------------
26! $Id: fft_xy.f90 1092 2013-02-02 11:24: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
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       INTEGER                             ::  sizw
600       COMPLEX, DIMENSION((nx+4)/2+1,nz+1) ::  work
601#endif
602
603       IF ( fft_method == 'temperton-algorithm' )  THEN
604
605          siza = SIZE( ai, 1 )
606
607          IF ( direction == 'forward')  THEN
608
609             ai(0:nx,1:nz) = ar(0:nx,1:nz)
610             ai(nx+1:,:)   = 0.0
611
612             CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
613
614             DO  k = 1, nz
615                DO  i = 0, (nx+1)/2
616                   ar(i,k) = ai(2*i,k)
617                ENDDO
618                DO  i = 1, (nx+1)/2 - 1
619                   ar(nx+1-i,k) = ai(2*i+1,k)
620                ENDDO
621             ENDDO
622
623          ELSE
624
625             DO  k = 1, nz
626                DO  i = 0, (nx+1)/2
627                   ai(2*i,k) = ar(i,k)
628                ENDDO
629                DO  i = 1, (nx+1)/2 - 1
630                   ai(2*i+1,k) = ar(nx+1-i,k)
631                ENDDO
632                ai(1,k) = 0.0
633                ai(nx+2,k) = 0.0
634             ENDDO
635
636             CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
637
638             ar(0:nx,1:nz) = ai(0:nx,1:nz)
639
640          ENDIF
641
642       ELSEIF ( fft_method == 'system-specific' )  THEN
643
644#if defined( __nec )
645          siza = SIZE( ai, 1 )
646          sizw = SIZE( work, 1 )
647
648          IF ( direction == 'forward')  THEN
649
650!
651!--          Tables are initialized once more. This call should not be
652!--          necessary, but otherwise program aborts in asymmetric case
653             CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
654                          trig_xf, work1, 0 )
655
656             ai(0:nx,1:nz) = ar(0:nx,1:nz)
657             IF ( nz1 > nz )  THEN
658                ai(:,nz1) = 0.0
659             ENDIF
660
661             CALL DZFFTM( 1, nx+1, nz1, sqr_nx, ai, siza, work, sizw, &
662                          trig_xf, work1, 0 )
663
664             DO  k = 1, nz
665                DO  i = 0, (nx+1)/2
666                   ar(i,k) = REAL( work(i+1,k) )
667                ENDDO
668                DO  i = 1, (nx+1)/2 - 1
669                   ar(nx+1-i,k) = AIMAG( work(i+1,k) )
670                ENDDO
671             ENDDO
672
673          ELSE
674
675!
676!--          Tables are initialized once more. This call should not be
677!--          necessary, but otherwise program aborts in asymmetric case
678             CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
679                          trig_xb, work1, 0 )
680
681             IF ( nz1 > nz )  THEN
682                work(:,nz1) = 0.0
683             ENDIF
684             DO  k = 1, nz
685                work(1,k) = CMPLX( ar(0,k), 0.0 )
686                DO  i = 1, (nx+1)/2 - 1
687                   work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
688                ENDDO
689                work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
690             ENDDO
691
692             CALL ZDFFTM( -1, nx+1, nz1, sqr_nx, work, sizw, ai, siza, &
693                          trig_xb, work1, 0 )
694
695             ar(0:nx,1:nz) = ai(0:nx,1:nz)
696
697          ENDIF
698
699#else
700          message_string = 'no system-specific fft-call available'
701          CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 ) 
702#endif
703
704       ELSE
705
706          message_string = 'fft method "' // TRIM( fft_method) // &
707                           '" not available'
708          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
709
710       ENDIF
711
712    END SUBROUTINE fft_x_m
713
714    SUBROUTINE fft_y_m( ar, ny1, direction )
715
716!----------------------------------------------------------------------!
717!                               fft_y_m                                !
718!                                                                      !
719!               Fourier-transformation along y-direction               !
720!                 Version for 1d domain decomposition                  !
721!             using multiple 1D FFT from Math Keisan on NEC            !
722!                       or Temperton-algorithm                         !
723!       (no singleton-algorithm on NEC because it does not vectorize)  !
724!                                                                      !
725!----------------------------------------------------------------------!
726
727       IMPLICIT NONE
728
729       CHARACTER (LEN=*)              ::  direction
730       INTEGER                        ::  j, k, ny1, siza
731
732       REAL, DIMENSION(0:ny1,nz)      ::  ar
733       REAL, DIMENSION(0:ny+3,nz+1)   ::  ai
734       REAL, DIMENSION(6*(ny+4),nz+1) ::  work1
735#if defined( __nec )
736       INTEGER                             ::  sizw
737       COMPLEX, DIMENSION((ny+4)/2+1,nz+1) ::  work
738#endif
739
740       IF ( fft_method == 'temperton-algorithm' )  THEN
741
742          siza = SIZE( ai, 1 )
743
744          IF ( direction == 'forward')  THEN
745
746             ai(0:ny,1:nz) = ar(0:ny,1:nz)
747             ai(ny+1:,:)   = 0.0
748
749             CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
750
751             DO  k = 1, nz
752                DO  j = 0, (ny+1)/2
753                   ar(j,k) = ai(2*j,k)
754                ENDDO
755                DO  j = 1, (ny+1)/2 - 1
756                   ar(ny+1-j,k) = ai(2*j+1,k)
757                ENDDO
758             ENDDO
759
760          ELSE
761
762             DO  k = 1, nz
763                DO  j = 0, (ny+1)/2
764                   ai(2*j,k) = ar(j,k)
765                ENDDO
766                DO  j = 1, (ny+1)/2 - 1
767                   ai(2*j+1,k) = ar(ny+1-j,k)
768                ENDDO
769                ai(1,k) = 0.0
770                ai(ny+2,k) = 0.0
771             ENDDO
772
773             CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
774
775             ar(0:ny,1:nz) = ai(0:ny,1:nz)
776
777          ENDIF
778
779       ELSEIF ( fft_method == 'system-specific' )  THEN
780
781#if defined( __nec )
782          siza = SIZE( ai, 1 )
783          sizw = SIZE( work, 1 )
784
785          IF ( direction == 'forward')  THEN
786
787!
788!--          Tables are initialized once more. This call should not be
789!--          necessary, but otherwise program aborts in asymmetric case
790             CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
791                          trig_yf, work1, 0 )
792
793             ai(0:ny,1:nz) = ar(0:ny,1:nz)
794             IF ( nz1 > nz )  THEN
795                ai(:,nz1) = 0.0
796             ENDIF
797
798             CALL DZFFTM( 1, ny+1, nz1, sqr_ny, ai, siza, work, sizw, &
799                          trig_yf, work1, 0 )
800
801             DO  k = 1, nz
802                DO  j = 0, (ny+1)/2
803                   ar(j,k) = REAL( work(j+1,k) )
804                ENDDO
805                DO  j = 1, (ny+1)/2 - 1
806                   ar(ny+1-j,k) = AIMAG( work(j+1,k) )
807                ENDDO
808             ENDDO
809
810          ELSE
811
812!
813!--          Tables are initialized once more. This call should not be
814!--          necessary, but otherwise program aborts in asymmetric case
815             CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
816                          trig_yb, work1, 0 )
817
818             IF ( nz1 > nz )  THEN
819                work(:,nz1) = 0.0
820             ENDIF
821             DO  k = 1, nz
822                work(1,k) = CMPLX( ar(0,k), 0.0 )
823                DO  j = 1, (ny+1)/2 - 1
824                   work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
825                ENDDO
826                work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
827             ENDDO
828
829             CALL ZDFFTM( -1, ny+1, nz1, sqr_ny, work, sizw, ai, siza, &
830                          trig_yb, work1, 0 )
831
832             ar(0:ny,1:nz) = ai(0:ny,1:nz)
833
834          ENDIF
835
836#else
837          message_string = 'no system-specific fft-call available'
838          CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 ) 
839#endif
840
841       ELSE
842         
843          message_string = 'fft method "' // TRIM( fft_method) // &
844                           '" not available'
845          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
846
847       ENDIF
848
849    END SUBROUTINE fft_y_m
850
851 END MODULE fft_xy
Note: See TracBrowser for help on using the repository browser.