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

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

code has been put under the GNU General Public License (v3)

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