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

Last change on this file since 1 was 1, checked in by raasch, 15 years ago

Initial repository layout and content

File size: 24.8 KB
Line 
1 MODULE fft_xy
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: fft_xy.f90,v $
11! Revision 1.4  2006/03/28 12:27:09  raasch
12! Stop when system-specific fft is selected on NEC. For unknown reasons this
13! causes a program abort during first allocation in init_grid.
14!
15! Revision 1.3  2006/02/23 12:22:41  raasch
16! System specfic routines are not allowed to be used on ibmy
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  2004/04/30 11:37:14  raasch
29! Initial revision
30!
31! Revision 1.6  2003/08/01 08:15:09  raasch
32! No abort in fft_1dd_init on t3e-systems if system specific routines are used
33!
34! Revision 1.5  2003/04/16 12:51:47  raasch
35! Temperton-fft implemented in routines fft_x_1dd and fft_y_1dd
36!
37! Revision 1.4  2003/03/16 09:37:26  raasch
38! Two underscores (_) are placed in front of all define-strings
39!
40! Revision 1.3  2003/03/12 16:29:39  raasch
41! Routines fft_x_1dd_m and fft_y_1dd_m added (suitable for multiple 1d-fft on
42! vector processors)
43!
44! Revision 1.2  2002/12/19 14:49:18  raasch
45! STOP statement replaced by call of subroutine local_stop
46!
47! Revision 1.1  2002/06/11 13:00:49  raasch
48! Initial revision
49!
50!
51! Description:
52! ------------
53! Fast Fourier transformation along x and y for 1d domain decomposition along x.
54! Original version: Klaus Ketelsen (May 2002)
55!------------------------------------------------------------------------------!
56
57    USE array_kind
58    USE control_parameters
59    USE indices
60    USE singleton
61    USE temperton_fft
62
63    IMPLICIT NONE
64
65    PRIVATE
66    PUBLIC fft_x, fft_y, fft_init, fft_x_m, fft_y_m
67
68    INTEGER, DIMENSION(:), ALLOCATABLE, SAVE ::  ifax_x, ifax_y
69
70    LOGICAL, SAVE                            ::  init_fft = .FALSE.
71
72    REAL, SAVE ::  sqr_nx, sqr_ny
73    REAL, DIMENSION(:), ALLOCATABLE, SAVE    ::  trigs_x, trigs_y
74
75#if defined( __ibm )
76    INTEGER, PARAMETER ::  nau1 = 20000, nau2 = 22000
77!
78!-- The following working arrays contain tables and have to be "save" and
79!-- shared in OpenMP sense
80    REAL, DIMENSION(nau1), SAVE ::  aux1, auy1, aux3, auy3
81#elif defined( __nec )
82    INTEGER, SAVE ::  nz1
83    REAL, DIMENSION(:), ALLOCATABLE, SAVE ::  trig_xb, trig_xf, trig_yb, &
84                                              trig_yf
85#endif
86
87!
88!-- Public interfaces
89    INTERFACE fft_init
90       MODULE PROCEDURE fft_init
91    END INTERFACE fft_init
92
93    INTERFACE fft_x
94       MODULE PROCEDURE fft_x
95    END INTERFACE fft_x
96
97    INTERFACE fft_y
98       MODULE PROCEDURE fft_y
99    END INTERFACE fft_y
100
101    INTERFACE fft_x_m
102       MODULE PROCEDURE fft_x_m
103    END INTERFACE fft_x_m
104
105    INTERFACE fft_y_m
106       MODULE PROCEDURE fft_y_m
107    END INTERFACE fft_y_m
108
109 CONTAINS
110
111
112    SUBROUTINE fft_init
113
114       IMPLICIT NONE
115
116!
117!--    The following temporary working arrays have to be on stack or private
118!--    in OpenMP sense
119#if defined( __ibm )
120       REAL, DIMENSION(0:nx+2) :: workx
121       REAL, DIMENSION(0:ny+2) :: worky
122       REAL, DIMENSION(nau2)   :: aux2, auy2, aux4, auy4
123#elif defined( __nec )
124       REAL, DIMENSION(0:nx+3,nz+1)   ::  work_x
125       REAL, DIMENSION(0:ny+3,nz+1)   ::  work_y
126       REAL, DIMENSION(6*(nx+3),nz+1) ::  workx
127       REAL, DIMENSION(6*(ny+3),nz+1) ::  worky
128#endif
129
130!
131!--    Return, if already called
132       IF ( init_fft )  THEN
133          RETURN
134       ELSE
135          init_fft = .TRUE.
136       ENDIF
137
138       IF ( fft_method == 'system-specific' )  THEN
139
140          sqr_nx = SQRT( 1.0 / ( nx + 1.0 ) )
141          sqr_ny = SQRT( 1.0 / ( ny + 1.0 ) )
142#if defined( __ibm )  &&  ! defined( __ibmy_special )
143!
144!--       Initialize tables for fft along x
145          CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1,  1, sqr_nx, aux1, nau1, &
146                      aux2, nau2 )
147          CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_nx, aux3, nau1, &
148                      aux4, nau2 )
149!
150!--       Initialize tables for fft along y
151          CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1,  1, sqr_ny, auy1, nau1, &
152                      auy2, nau2 )
153          CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_ny, auy3, nau1, &
154                      auy4, nau2 )
155#elif defined( __nec )
156          PRINT*, '+++ fft_init: fft method "', fft_method, &
157                  '" currently does not work on NEC'
158          CALL local_stop
159
160          ALLOCATE( trig_xb(2*(nx+1)), trig_xf(2*(nx+1)), &
161                    trig_yb(2*(ny+1)), trig_yf(2*(ny+1)) )
162
163          work_x = 0.0
164          work_y = 0.0
165          nz1  = nz + MOD( nz+1, 2 )  ! odd nz slows down fft significantly
166                                      ! when using the NEC ffts
167
168!
169!--       Initialize tables for fft along x (non-vector and vector case (M))
170          CALL DZFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xf, workx, 0 )
171          CALL ZDFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xb, workx, 0 )
172          CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
173                       trig_xf, workx, 0 )
174          CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
175                       trig_xb, workx, 0 )
176!
177!--       Initialize tables for fft along y (non-vector and vector case (M))
178          CALL DZFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yf, worky, 0 )
179          CALL ZDFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yb, worky, 0 )
180          CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
181                       trig_yf, worky, 0 )
182          CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
183                       trig_yb, worky, 0 )
184#else
185          PRINT*, '+++ fft_init: no system-specific fft-call available'
186          CALL local_stop
187#endif
188       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
189!
190!--       Temperton-algorithm
191!--       Initialize tables for fft along x and y
192          ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
193
194          CALL set99( trigs_x, ifax_x, nx+1 )
195          CALL set99( trigs_y, ifax_y, ny+1 )
196
197       ELSEIF ( fft_method == 'singleton-algorithm' )  THEN
198
199          CONTINUE
200
201       ELSE
202
203          PRINT*, '+++ fft_init: fft method "', fft_method, &
204                  '" not available'
205          CALL local_stop
206
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          PRINT*, '+++ fft_x: no system-specific fft-call available'
380          CALL local_stop
381#endif
382       ELSE
383
384          PRINT*, '+++ fft_x: fft method "', fft_method, '" not available'
385          CALL local_stop
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          PRINT*, '+++ fft_y: no system-specific fft-call available'
560          CALL local_stop
561#endif
562
563       ELSE
564
565          PRINT*, '+++ fft_y: fft method "', fft_method, '" not available'
566          CALL local_stop
567
568       ENDIF
569
570    END SUBROUTINE fft_y
571
572    SUBROUTINE fft_x_m( ar, direction )
573
574!----------------------------------------------------------------------!
575!                               fft_x_m                                !
576!                                                                      !
577!               Fourier-transformation along x-direction               !
578!                 Version for 1d domain decomposition                  !
579!             using multiple 1D FFT from Math Keisan on NEC            !
580!                       or Temperton-algorithm                         !
581!       (no singleton-algorithm on NEC because it does not vectorize)  !
582!                                                                      !
583!----------------------------------------------------------------------!
584
585       IMPLICIT NONE
586
587       CHARACTER (LEN=*)              ::  direction
588       INTEGER                        ::  i, k, siza, sizw
589
590       REAL, DIMENSION(0:nx,nz)       ::  ar
591       REAL, DIMENSION(0:nx+3,nz+1)   ::  ai
592       REAL, DIMENSION(6*(nx+4),nz+1) ::  work1
593#if defined( __nec )
594       COMPLEX, DIMENSION((nx+4)/2+1,nz+1) ::  work
595#endif
596
597       IF ( fft_method == 'temperton-algorithm' )  THEN
598
599          siza = SIZE( ai, 1 )
600
601          IF ( direction == 'forward')  THEN
602
603             ai(0:nx,1:nz) = ar(0:nx,1:nz)
604             ai(nx+1:,:)   = 0.0
605
606             CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
607
608             DO  k = 1, nz
609                DO  i = 0, (nx+1)/2
610                   ar(i,k) = ai(2*i,k)
611                ENDDO
612                DO  i = 1, (nx+1)/2 - 1
613                   ar(nx+1-i,k) = ai(2*i+1,k)
614                ENDDO
615             ENDDO
616
617          ELSE
618
619             DO  k = 1, nz
620                DO  i = 0, (nx+1)/2
621                   ai(2*i,k) = ar(i,k)
622                ENDDO
623                DO  i = 1, (nx+1)/2 - 1
624                   ai(2*i+1,k) = ar(nx+1-i,k)
625                ENDDO
626                ai(1,k) = 0.0
627                ai(nx+2,k) = 0.0
628             ENDDO
629
630             CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
631
632             ar(0:nx,1:nz) = ai(0:nx,1:nz)
633
634          ENDIF
635
636       ELSEIF ( fft_method == 'system-specific' )  THEN
637
638#if defined( __nec )
639          siza = SIZE( ai, 1 )
640          sizw = SIZE( work, 1 )
641
642          IF ( direction == 'forward')  THEN
643
644!
645!--          Tables are initialized once more. This call should not be
646!--          necessary, but otherwise program aborts in asymmetric case
647             CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
648                          trig_xf, work1, 0 )
649
650             ai(0:nx,1:nz) = ar(0:nx,1:nz)
651             IF ( nz1 > nz )  THEN
652                ai(:,nz1) = 0.0
653             ENDIF
654
655             CALL DZFFTM( 1, nx+1, nz1, sqr_nx, ai, siza, work, sizw, &
656                          trig_xf, work1, 0 )
657
658             DO  k = 1, nz
659                DO  i = 0, (nx+1)/2
660                   ar(i,k) = REAL( work(i+1,k) )
661                ENDDO
662                DO  i = 1, (nx+1)/2 - 1
663                   ar(nx+1-i,k) = AIMAG( work(i+1,k) )
664                ENDDO
665             ENDDO
666
667          ELSE
668
669!
670!--          Tables are initialized once more. This call should not be
671!--          necessary, but otherwise program aborts in asymmetric case
672             CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
673                          trig_xb, work1, 0 )
674
675             IF ( nz1 > nz )  THEN
676                work(:,nz1) = 0.0
677             ENDIF
678             DO  k = 1, nz
679                work(1,k) = CMPLX( ar(0,k), 0.0 )
680                DO  i = 1, (nx+1)/2 - 1
681                   work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
682                ENDDO
683                work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
684             ENDDO
685
686             CALL ZDFFTM( -1, nx+1, nz1, sqr_nx, work, sizw, ai, siza, &
687                          trig_xb, work1, 0 )
688
689             ar(0:nx,1:nz) = ai(0:nx,1:nz)
690
691          ENDIF
692
693#else
694          PRINT*, '+++ fft_x_m: no system-specific fft-call available'
695          STOP
696#endif
697
698       ELSE
699
700          PRINT*, '+++ fft_x_m: fft method "', fft_method, '" not available'
701          CALL local_stop
702
703       ENDIF
704
705    END SUBROUTINE fft_x_m
706
707    SUBROUTINE fft_y_m( ar, ny1, direction )
708
709!----------------------------------------------------------------------!
710!                               fft_y_m                                !
711!                                                                      !
712!               Fourier-transformation along y-direction               !
713!                 Version for 1d domain decomposition                  !
714!             using multiple 1D FFT from Math Keisan on NEC            !
715!                       or Temperton-algorithm                         !
716!       (no singleton-algorithm on NEC because it does not vectorize)  !
717!                                                                      !
718!----------------------------------------------------------------------!
719
720       IMPLICIT NONE
721
722       CHARACTER (LEN=*)              ::  direction
723       INTEGER                        ::  j, k, ny1, siza, sizw
724
725       REAL, DIMENSION(0:ny1,nz)      ::  ar
726       REAL, DIMENSION(0:ny+3,nz+1)   ::  ai
727       REAL, DIMENSION(6*(ny+4),nz+1) ::  work1
728#if defined( __nec )
729       COMPLEX, DIMENSION((ny+4)/2+1,nz+1) ::  work
730#endif
731
732       IF ( fft_method == 'temperton-algorithm' )  THEN
733
734          siza = SIZE( ai, 1 )
735
736          IF ( direction == 'forward')  THEN
737
738             ai(0:ny,1:nz) = ar(0:ny,1:nz)
739             ai(ny+1:,:)   = 0.0
740
741             CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
742
743             DO  k = 1, nz
744                DO  j = 0, (ny+1)/2
745                   ar(j,k) = ai(2*j,k)
746                ENDDO
747                DO  j = 1, (ny+1)/2 - 1
748                   ar(ny+1-j,k) = ai(2*j+1,k)
749                ENDDO
750             ENDDO
751
752          ELSE
753
754             DO  k = 1, nz
755                DO  j = 0, (ny+1)/2
756                   ai(2*j,k) = ar(j,k)
757                ENDDO
758                DO  j = 1, (ny+1)/2 - 1
759                   ai(2*j+1,k) = ar(ny+1-j,k)
760                ENDDO
761                ai(1,k) = 0.0
762                ai(ny+2,k) = 0.0
763             ENDDO
764
765             CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
766
767             ar(0:ny,1:nz) = ai(0:ny,1:nz)
768
769          ENDIF
770
771       ELSEIF ( fft_method == 'system-specific' )  THEN
772
773#if defined( __nec )
774          siza = SIZE( ai, 1 )
775          sizw = SIZE( work, 1 )
776
777          IF ( direction == 'forward')  THEN
778
779!
780!--          Tables are initialized once more. This call should not be
781!--          necessary, but otherwise program aborts in asymmetric case
782             CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
783                          trig_yf, work1, 0 )
784
785             ai(0:ny,1:nz) = ar(0:ny,1:nz)
786             IF ( nz1 > nz )  THEN
787                ai(:,nz1) = 0.0
788             ENDIF
789
790             CALL DZFFTM( 1, ny+1, nz1, sqr_ny, ai, siza, work, sizw, &
791                          trig_yf, work1, 0 )
792
793             DO  k = 1, nz
794                DO  j = 0, (ny+1)/2
795                   ar(j,k) = REAL( work(j+1,k) )
796                ENDDO
797                DO  j = 1, (ny+1)/2 - 1
798                   ar(ny+1-j,k) = AIMAG( work(j+1,k) )
799                ENDDO
800             ENDDO
801
802          ELSE
803
804!
805!--          Tables are initialized once more. This call should not be
806!--          necessary, but otherwise program aborts in asymmetric case
807             CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
808                          trig_yb, work1, 0 )
809
810             IF ( nz1 > nz )  THEN
811                work(:,nz1) = 0.0
812             ENDIF
813             DO  k = 1, nz
814                work(1,k) = CMPLX( ar(0,k), 0.0 )
815                DO  j = 1, (ny+1)/2 - 1
816                   work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
817                ENDDO
818                work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
819             ENDDO
820
821             CALL ZDFFTM( -1, ny+1, nz1, sqr_ny, work, sizw, ai, siza, &
822                          trig_yb, work1, 0 )
823
824             ar(0:ny,1:nz) = ai(0:ny,1:nz)
825
826          ENDIF
827
828#else
829          PRINT*, '+++ fft_y_m: no system-specific fft-call available'
830          STOP
831#endif
832
833       ELSE
834
835          PRINT*, '+++ fft_y_m: fft method "', fft_method, '" not available'
836          CALL local_stop
837
838       ENDIF
839
840    END SUBROUTINE fft_y_m
841
842 END MODULE fft_xy
Note: See TracBrowser for help on using the repository browser.