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

Last change on this file since 274 was 274, checked in by heinze, 12 years ago

Indentation of the message calls corrected

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