source: palm/trunk/SOURCE/fft_xy_mod.f90 @ 2118

Last change on this file since 2118 was 2118, checked in by raasch, 4 years ago

all OpenACC directives and related parts removed from the code

  • Property svn:keywords set to Id
File size: 52.9 KB
Line 
1!> @file fft_xy_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! 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-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22! OpenACC directives and CUDA-fft related code removed
23!
24! Former revisions:
25! -----------------
26! $Id: fft_xy_mod.f90 2118 2017-01-17 16:38:49Z raasch $
27!
28! 2000 2016-08-20 18:09:15Z knoop
29! Forced header and separation lines into 80 columns
30!
31! 1850 2016-04-08 13:29:27Z maronga
32! Module renamed
33!
34! 1815 2016-04-06 13:49:59Z raasch
35! cpp-directives for ibmy removed
36!
37! 1749 2016-02-09 12:19:56Z raasch
38! small OpenACC bugfix
39!
40! 1682 2015-10-07 23:56:08Z knoop
41! Code annotations made doxygen readable
42!
43! 1600 2015-06-11 15:50:12Z raasch
44! bugfix: openMP threadprivate statement moved after variable declaration
45!
46! 1482 2014-10-18 12:34:45Z raasch
47! cudafft workaround for data declaration of ar_tmp because of PGI 14.1 bug
48!
49! 1402 2014-05-09 14:25:13Z raasch
50! fortran bugfix for r1392
51!
52! 1398 2014-05-07 11:15:00Z heinze
53! bugfix: typo removed for KIND in CMPLX function
54!
55! 1392 2014-05-06 09:10:05Z raasch
56! bugfix: KIND attribute added to CMPLX functions
57!
58! 1374 2014-04-25 12:55:07Z raasch
59! bugfixes: missing variables added to ONLY list, dpk renamed dp
60!
61! 1372 2014-04-24 06:29:32Z raasch
62! openMP-bugfix for fftw: some arrays defined as threadprivate
63!
64! 1353 2014-04-08 15:21:23Z heinze
65! REAL constants provided with KIND-attribute
66!
67! 1342 2014-03-26 17:04:47Z kanani
68! REAL constants defined as wp-kind
69!
70! 1322 2014-03-20 16:38:49Z raasch
71! REAL functions provided with KIND-attribute
72!
73! 1320 2014-03-20 08:40:49Z raasch
74! ONLY-attribute added to USE-statements,
75! kind-parameters added to all INTEGER and REAL declaration statements,
76! kinds are defined in new module kinds,
77! old module precision_kind is removed,
78! revision history before 2012 removed,
79! comment fields (!:) to be used for variable explanations added to
80! all variable declaration statements
81!
82! 1304 2014-03-12 10:29:42Z raasch
83! openmp bugfix: work1 used in Temperton algorithm must be private
84!
85! 1257 2013-11-08 15:18:40Z raasch
86! openacc loop and loop vector clauses removed, declare create moved after
87! the FORTRAN declaration statement
88!
89! 1219 2013-08-30 09:33:18Z heinze
90! bugfix: use own branch for fftw
91!
92! 1216 2013-08-26 09:31:42Z raasch
93! fft_x and fft_y modified for parallel / ovverlapping execution of fft and
94! transpositions,
95! fftw implemented for 1d-decomposition (fft_x_1d, fft_y_1d)
96!
97! 1210 2013-08-14 10:58:20Z raasch
98! fftw added
99!
100! 1166 2013-05-24 13:55:44Z raasch
101! C_DOUBLE/COMPLEX reset to dpk
102!
103! 1153 2013-05-10 14:33:08Z raasch
104! code adjustment of data types for CUDA fft required by PGI 12.3 / CUDA 5.0
105!
106! 1111 2013-03-08 23:54:10Z raasch
107! further openACC statements added, CUDA branch completely runs on GPU
108! bugfix: CUDA fft plans adjusted for domain decomposition (before they always
109! used total domain)
110!
111! 1106 2013-03-04 05:31:38Z raasch
112! CUDA fft added
113! array_kind renamed precision_kind, 3D- instead of 1D-loops in fft_x and fft_y
114! old fft_x, fft_y become fft_x_1d, fft_y_1d and are used for 1D-decomposition
115!
116! 1092 2013-02-02 11:24:22Z raasch
117! variable sizw declared for NEC case only
118!
119! 1036 2012-10-22 13:43:42Z raasch
120! code put under GPL (PALM 3.9)
121!
122! Revision 1.1  2002/06/11 13:00:49  raasch
123! Initial revision
124!
125!
126! Description:
127! ------------
128!> Fast Fourier transformation along x and y for 1d domain decomposition along x.
129!> Original version: Klaus Ketelsen (May 2002)
130!------------------------------------------------------------------------------!
131 MODULE fft_xy
132 
133
134    USE control_parameters,                                                    &
135        ONLY:  fft_method, message_string
136       
137    USE indices,                                                               &
138        ONLY:  nx, ny, nz
139       
140#if defined( __fftw )
141    USE, INTRINSIC ::  ISO_C_BINDING
142#endif
143
144    USE kinds
145   
146    USE singleton,                                                             &
147        ONLY: fftn
148   
149    USE temperton_fft
150   
151    USE transpose_indices,                                                     &
152        ONLY:  nxl_y, nxr_y, nyn_x, nys_x, nzb_x, nzb_y, nzt_x, nzt_y
153
154    IMPLICIT NONE
155
156    PRIVATE
157    PUBLIC fft_x, fft_x_1d, fft_y, fft_y_1d, fft_init, fft_x_m, fft_y_m
158
159    INTEGER(iwp), DIMENSION(:), ALLOCATABLE, SAVE ::  ifax_x  !<
160    INTEGER(iwp), DIMENSION(:), ALLOCATABLE, SAVE ::  ifax_y  !<
161
162    LOGICAL, SAVE ::  init_fft = .FALSE.  !<
163
164    REAL(wp), SAVE ::  dnx      !<
165    REAL(wp), SAVE ::  dny      !<
166    REAL(wp), SAVE ::  sqr_dnx  !<
167    REAL(wp), SAVE ::  sqr_dny  !<
168   
169    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trigs_x  !<
170    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trigs_y  !<
171
172#if defined( __ibm )
173    INTEGER(iwp), PARAMETER ::  nau1 = 20000  !<
174    INTEGER(iwp), PARAMETER ::  nau2 = 22000  !<
175!
176!-- The following working arrays contain tables and have to be "save" and
177!-- shared in OpenMP sense
178    REAL(wp), DIMENSION(nau1), SAVE ::  aux1  !<
179    REAL(wp), DIMENSION(nau1), SAVE ::  auy1  !<
180    REAL(wp), DIMENSION(nau1), SAVE ::  aux3  !<
181    REAL(wp), DIMENSION(nau1), SAVE ::  auy3  !<
182   
183#elif defined( __nec )
184    INTEGER(iwp), SAVE ::  nz1  !<
185   
186    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trig_xb  !<
187    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trig_xf  !<
188    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trig_yb  !<
189    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trig_yf  !<
190   
191#endif
192
193#if defined( __fftw )
194    INCLUDE  'fftw3.f03'
195    INTEGER(KIND=C_INT) ::  nx_c  !<
196    INTEGER(KIND=C_INT) ::  ny_c  !<
197   
198    COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE ::  x_out  !<
199    COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE ::         &
200       y_out  !<
201   
202    REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE ::                    &
203       x_in   !<
204    REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE ::                    &
205       y_in   !<
206    !$OMP THREADPRIVATE( x_out, y_out, x_in, y_in )
207   
208   
209    TYPE(C_PTR), SAVE ::  plan_xf, plan_xi, plan_yf, plan_yi
210#endif
211
212!
213!-- Public interfaces
214    INTERFACE fft_init
215       MODULE PROCEDURE fft_init
216    END INTERFACE fft_init
217
218    INTERFACE fft_x
219       MODULE PROCEDURE fft_x
220    END INTERFACE fft_x
221
222    INTERFACE fft_x_1d
223       MODULE PROCEDURE fft_x_1d
224    END INTERFACE fft_x_1d
225
226    INTERFACE fft_y
227       MODULE PROCEDURE fft_y
228    END INTERFACE fft_y
229
230    INTERFACE fft_y_1d
231       MODULE PROCEDURE fft_y_1d
232    END INTERFACE fft_y_1d
233
234    INTERFACE fft_x_m
235       MODULE PROCEDURE fft_x_m
236    END INTERFACE fft_x_m
237
238    INTERFACE fft_y_m
239       MODULE PROCEDURE fft_y_m
240    END INTERFACE fft_y_m
241
242 CONTAINS
243
244
245!------------------------------------------------------------------------------!
246! Description:
247! ------------
248!> @todo Missing subroutine description.
249!------------------------------------------------------------------------------!
250    SUBROUTINE fft_init
251
252       IMPLICIT NONE
253
254!
255!--    The following temporary working arrays have to be on stack or private
256!--    in OpenMP sense
257#if defined( __ibm )
258       REAL(wp), DIMENSION(0:nx+2) ::  workx  !<
259       REAL(wp), DIMENSION(0:ny+2) ::  worky  !<
260       REAL(wp), DIMENSION(nau2)   ::  aux2   !<
261       REAL(wp), DIMENSION(nau2)   ::  auy2   !<
262       REAL(wp), DIMENSION(nau2)   ::  aux4   !<
263       REAL(wp), DIMENSION(nau2)   ::  auy4   !<
264#elif defined( __nec )
265       REAL(wp), DIMENSION(0:nx+3,nz+1)   ::  work_x  !<
266       REAL(wp), DIMENSION(0:ny+3,nz+1)   ::  work_y  !<
267       REAL(wp), DIMENSION(6*(nx+3),nz+1) ::  workx   !<
268       REAL(wp), DIMENSION(6*(ny+3),nz+1) ::  worky   !<
269#endif
270
271!
272!--    Return, if already called
273       IF ( init_fft )  THEN
274          RETURN
275       ELSE
276          init_fft = .TRUE.
277       ENDIF
278
279       IF ( fft_method == 'system-specific' )  THEN
280
281          dnx = 1.0_wp / ( nx + 1.0_wp )
282          dny = 1.0_wp / ( ny + 1.0_wp )
283          sqr_dnx = SQRT( dnx )
284          sqr_dny = SQRT( dny )
285#if defined( __ibm )
286!
287!--       Initialize tables for fft along x
288          CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1,  1, sqr_dnx, aux1, nau1, &
289                      aux2, nau2 )
290          CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
291                      aux4, nau2 )
292!
293!--       Initialize tables for fft along y
294          CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1,  1, sqr_dny, auy1, nau1, &
295                      auy2, nau2 )
296          CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
297                      auy4, nau2 )
298#elif defined( __nec )
299          message_string = 'fft method "' // TRIM( fft_method) // &
300                           '" currently does not work on NEC'
301          CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 )
302
303          ALLOCATE( trig_xb(2*(nx+1)), trig_xf(2*(nx+1)),                      &
304                    trig_yb(2*(ny+1)), trig_yf(2*(ny+1)) )
305
306          work_x = 0.0_wp
307          work_y = 0.0_wp
308          nz1  = nz + MOD( nz+1, 2 )  ! odd nz slows down fft significantly
309                                      ! when using the NEC ffts
310
311!
312!--       Initialize tables for fft along x (non-vector and vector case (M))
313          CALL DZFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xf, workx, 0 )
314          CALL ZDFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xb, workx, 0 )
315          CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4,      &
316                       trig_xf, workx, 0 )
317          CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4,      &
318                       trig_xb, workx, 0 )
319!
320!--       Initialize tables for fft along y (non-vector and vector case (M))
321          CALL DZFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yf, worky, 0 )
322          CALL ZDFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yb, worky, 0 )
323          CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4,      &
324                       trig_yf, worky, 0 )
325          CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4,      &
326                       trig_yb, worky, 0 )
327#else
328          message_string = 'no system-specific fft-call available'
329          CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
330#endif
331       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
332!
333!--       Temperton-algorithm
334!--       Initialize tables for fft along x and y
335          ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
336
337          CALL set99( trigs_x, ifax_x, nx+1 )
338          CALL set99( trigs_y, ifax_y, ny+1 )
339
340       ELSEIF ( fft_method == 'fftw' )  THEN
341!
342!--       FFTW
343#if defined( __fftw )
344          nx_c = nx+1
345          ny_c = ny+1
346          !$OMP PARALLEL
347          ALLOCATE( x_in(0:nx+2), y_in(0:ny+2), x_out(0:(nx+1)/2),             &
348                    y_out(0:(ny+1)/2) )
349          !$OMP END PARALLEL
350          plan_xf = FFTW_PLAN_DFT_R2C_1D( nx_c, x_in, x_out, FFTW_ESTIMATE )
351          plan_xi = FFTW_PLAN_DFT_C2R_1D( nx_c, x_out, x_in, FFTW_ESTIMATE )
352          plan_yf = FFTW_PLAN_DFT_R2C_1D( ny_c, y_in, y_out, FFTW_ESTIMATE )
353          plan_yi = FFTW_PLAN_DFT_C2R_1D( ny_c, y_out, y_in, FFTW_ESTIMATE )
354#else
355          message_string = 'preprocessor switch for fftw is missing'
356          CALL message( 'fft_init', 'PA0080', 1, 2, 0, 6, 0 )
357#endif
358
359       ELSEIF ( fft_method == 'singleton-algorithm' )  THEN
360
361          CONTINUE
362
363       ELSE
364
365          message_string = 'fft method "' // TRIM( fft_method) // &
366                           '" not available'
367          CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
368       ENDIF
369
370    END SUBROUTINE fft_init
371
372
373!------------------------------------------------------------------------------!
374! Description:
375! ------------
376!> Fourier-transformation along x-direction.                 
377!> Version for 2D-decomposition.
378!> It uses internal algorithms (Singleton or Temperton) or     
379!> system-specific routines, if they are available           
380!------------------------------------------------------------------------------!
381 
382    SUBROUTINE fft_x( ar, direction, ar_2d )
383
384
385       IMPLICIT NONE
386
387       CHARACTER (LEN=*) ::  direction  !<
388       
389       COMPLEX(wp), DIMENSION(:), ALLOCATABLE ::  cwork  !<
390
391       INTEGER(iwp) ::  i          !<
392       INTEGER(iwp) ::  ishape(1)  !<
393       INTEGER(iwp) ::  j          !<
394       INTEGER(iwp) ::  k          !<
395
396       LOGICAL ::  forward_fft !<
397       
398       REAL(wp), DIMENSION(0:nx+2) ::  work   !<
399       REAL(wp), DIMENSION(nx+2)   ::  work1  !<
400       
401#if defined( __ibm )
402       REAL(wp), DIMENSION(nau2) ::  aux2  !<
403       REAL(wp), DIMENSION(nau2) ::  aux4  !<
404#elif defined( __nec )
405       REAL(wp), DIMENSION(6*(nx+1)) ::  work2  !<
406#endif
407
408       REAL(wp), DIMENSION(0:nx,nys_x:nyn_x), OPTIONAL   ::                    &
409          ar_2d   !<
410       REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::                    &
411          ar      !<
412
413       IF ( direction == 'forward' )  THEN
414          forward_fft = .TRUE.
415       ELSE
416          forward_fft = .FALSE.
417       ENDIF
418
419       IF ( fft_method == 'singleton-algorithm' )  THEN
420
421!
422!--       Performing the fft with singleton's software works on every system,
423!--       since it is part of the model
424          ALLOCATE( cwork(0:nx) )
425     
426          IF ( forward_fft )   then
427
428             !$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
429             !$OMP DO
430             DO  k = nzb_x, nzt_x
431                DO  j = nys_x, nyn_x
432
433                   DO  i = 0, nx
434                      cwork(i) = CMPLX( ar(i,j,k), KIND=wp )
435                   ENDDO
436
437                   ishape = SHAPE( cwork )
438                   CALL FFTN( cwork, ishape )
439
440                   DO  i = 0, (nx+1)/2
441                      ar(i,j,k) = REAL( cwork(i), KIND=wp )
442                   ENDDO
443                   DO  i = 1, (nx+1)/2 - 1
444                      ar(nx+1-i,j,k) = -AIMAG( cwork(i) )
445                   ENDDO
446
447                ENDDO
448             ENDDO
449             !$OMP END PARALLEL
450
451          ELSE
452
453             !$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
454             !$OMP DO
455             DO  k = nzb_x, nzt_x
456                DO  j = nys_x, nyn_x
457
458                   cwork(0) = CMPLX( ar(0,j,k), 0.0_wp, KIND=wp )
459                   DO  i = 1, (nx+1)/2 - 1
460                      cwork(i)      = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k),       &
461                                             KIND=wp )
462                      cwork(nx+1-i) = CMPLX( ar(i,j,k),  ar(nx+1-i,j,k),       &
463                                             KIND=wp )
464                   ENDDO
465                   cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp, KIND=wp )
466
467                   ishape = SHAPE( cwork )
468                   CALL FFTN( cwork, ishape, inv = .TRUE. )
469
470                   DO  i = 0, nx
471                      ar(i,j,k) = REAL( cwork(i), KIND=wp )
472                   ENDDO
473
474                ENDDO
475             ENDDO
476             !$OMP END PARALLEL
477
478          ENDIF
479
480          DEALLOCATE( cwork )
481
482       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
483
484!
485!--       Performing the fft with Temperton's software works on every system,
486!--       since it is part of the model
487          IF ( forward_fft )  THEN
488
489             !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
490             !$OMP DO
491             DO  k = nzb_x, nzt_x
492                DO  j = nys_x, nyn_x
493
494                   work(0:nx) = ar(0:nx,j,k)
495                   CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
496
497                   DO  i = 0, (nx+1)/2
498                      ar(i,j,k) = work(2*i)
499                   ENDDO
500                   DO  i = 1, (nx+1)/2 - 1
501                      ar(nx+1-i,j,k) = work(2*i+1)
502                   ENDDO
503
504                ENDDO
505             ENDDO
506             !$OMP END PARALLEL
507
508          ELSE
509
510             !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
511             !$OMP DO
512             DO  k = nzb_x, nzt_x
513                DO  j = nys_x, nyn_x
514
515                   DO  i = 0, (nx+1)/2
516                      work(2*i) = ar(i,j,k)
517                   ENDDO
518                   DO  i = 1, (nx+1)/2 - 1
519                      work(2*i+1) = ar(nx+1-i,j,k)
520                   ENDDO
521                   work(1)    = 0.0_wp
522                   work(nx+2) = 0.0_wp
523
524                   CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
525                   ar(0:nx,j,k) = work(0:nx)
526
527                ENDDO
528             ENDDO
529             !$OMP END PARALLEL
530
531          ENDIF
532
533       ELSEIF ( fft_method == 'fftw' )  THEN
534
535#if defined( __fftw )
536          IF ( forward_fft )  THEN
537
538             !$OMP PARALLEL PRIVATE ( work, i, j, k )
539             !$OMP DO
540             DO  k = nzb_x, nzt_x
541                DO  j = nys_x, nyn_x
542
543                   x_in(0:nx) = ar(0:nx,j,k)
544                   CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
545
546                   IF ( PRESENT( ar_2d ) )  THEN
547
548                      DO  i = 0, (nx+1)/2
549                         ar_2d(i,j) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
550                      ENDDO
551                      DO  i = 1, (nx+1)/2 - 1
552                         ar_2d(nx+1-i,j) = AIMAG( x_out(i) ) / ( nx+1 )
553                      ENDDO
554
555                   ELSE
556
557                      DO  i = 0, (nx+1)/2
558                         ar(i,j,k) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
559                      ENDDO
560                      DO  i = 1, (nx+1)/2 - 1
561                         ar(nx+1-i,j,k) = AIMAG( x_out(i) ) / ( nx+1 )
562                      ENDDO
563
564                   ENDIF
565
566                ENDDO
567             ENDDO
568             !$OMP END PARALLEL
569
570          ELSE
571             !$OMP PARALLEL PRIVATE ( work, i, j, k )
572             !$OMP DO
573             DO  k = nzb_x, nzt_x
574                DO  j = nys_x, nyn_x
575
576                   IF ( PRESENT( ar_2d ) )  THEN
577
578                      x_out(0) = CMPLX( ar_2d(0,j), 0.0_wp, KIND=wp )
579                      DO  i = 1, (nx+1)/2 - 1
580                         x_out(i) = CMPLX( ar_2d(i,j), ar_2d(nx+1-i,j),        &
581                                           KIND=wp )
582                      ENDDO
583                      x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0_wp,      &
584                                               KIND=wp )
585
586                   ELSE
587
588                      x_out(0) = CMPLX( ar(0,j,k), 0.0_wp, KIND=wp )
589                      DO  i = 1, (nx+1)/2 - 1
590                         x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k), KIND=wp )
591                      ENDDO
592                      x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp,       &
593                                               KIND=wp )
594
595                   ENDIF
596
597                   CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
598                   ar(0:nx,j,k) = x_in(0:nx)
599
600                ENDDO
601             ENDDO
602             !$OMP END PARALLEL
603
604          ENDIF
605#endif
606
607       ELSEIF ( fft_method == 'system-specific' )  THEN
608
609#if defined( __ibm )
610          IF ( forward_fft )  THEN
611
612             !$OMP PARALLEL PRIVATE ( work, i, j, k )
613             !$OMP DO
614             DO  k = nzb_x, nzt_x
615                DO  j = nys_x, nyn_x
616
617                   CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1,   &
618                               nau1, aux2, nau2 )
619
620                   DO  i = 0, (nx+1)/2
621                      ar(i,j,k) = work(2*i)
622                   ENDDO
623                   DO  i = 1, (nx+1)/2 - 1
624                      ar(nx+1-i,j,k) = work(2*i+1)
625                   ENDDO
626
627                ENDDO
628             ENDDO
629             !$OMP END PARALLEL
630
631          ELSE
632
633             !$OMP PARALLEL PRIVATE ( work, i, j, k )
634             !$OMP DO
635             DO  k = nzb_x, nzt_x
636                DO  j = nys_x, nyn_x
637
638                   DO  i = 0, (nx+1)/2
639                      work(2*i) = ar(i,j,k)
640                   ENDDO
641                   DO  i = 1, (nx+1)/2 - 1
642                      work(2*i+1) = ar(nx+1-i,j,k)
643                   ENDDO
644                   work(1) = 0.0_wp
645                   work(nx+2) = 0.0_wp
646
647                   CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx,      & 
648                               aux3, nau1, aux4, nau2 )
649
650                   DO  i = 0, nx
651                      ar(i,j,k) = work(i)
652                   ENDDO
653
654                ENDDO
655             ENDDO
656             !$OMP END PARALLEL
657
658          ENDIF
659
660#elif defined( __nec )
661
662          IF ( forward_fft )  THEN
663
664             !$OMP PARALLEL PRIVATE ( work, i, j, k )
665             !$OMP DO
666             DO  k = nzb_x, nzt_x
667                DO  j = nys_x, nyn_x
668
669                   work(0:nx) = ar(0:nx,j,k)
670
671                   CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
672     
673                   DO  i = 0, (nx+1)/2
674                      ar(i,j,k) = work(2*i)
675                   ENDDO
676                   DO  i = 1, (nx+1)/2 - 1
677                      ar(nx+1-i,j,k) = work(2*i+1)
678                   ENDDO
679
680                ENDDO
681             ENDDO
682             !$END OMP PARALLEL
683
684          ELSE
685
686             !$OMP PARALLEL PRIVATE ( work, i, j, k )
687             !$OMP DO
688             DO  k = nzb_x, nzt_x
689                DO  j = nys_x, nyn_x
690
691                   DO  i = 0, (nx+1)/2
692                      work(2*i) = ar(i,j,k)
693                   ENDDO
694                   DO  i = 1, (nx+1)/2 - 1
695                      work(2*i+1) = ar(nx+1-i,j,k)
696                   ENDDO
697                   work(1) = 0.0_wp
698                   work(nx+2) = 0.0_wp
699
700                   CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
701
702                   ar(0:nx,j,k) = work(0:nx)
703
704                ENDDO
705             ENDDO
706             !$OMP END PARALLEL
707
708          ENDIF
709
710#else
711          message_string = 'no system-specific fft-call available'
712          CALL message( 'fft_x', 'PA0188', 1, 2, 0, 6, 0 )
713#endif
714
715       ELSE
716
717          message_string = 'fft method "' // TRIM( fft_method) // &
718                           '" not available'
719          CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
720
721       ENDIF
722
723    END SUBROUTINE fft_x
724
725!------------------------------------------------------------------------------!
726! Description:
727! ------------
728!> Fourier-transformation along x-direction.
729!> Version for 1D-decomposition.
730!> It uses internal algorithms (Singleton or Temperton) or
731!> system-specific routines, if they are available
732!------------------------------------------------------------------------------!
733 
734    SUBROUTINE fft_x_1d( ar, direction )
735
736
737       IMPLICIT NONE
738
739       CHARACTER (LEN=*) ::  direction  !<
740       
741       INTEGER(iwp) ::  i               !<
742       INTEGER(iwp) ::  ishape(1)       !<
743
744       LOGICAL ::  forward_fft          !<
745
746       REAL(wp), DIMENSION(0:nx)   ::  ar     !<
747       REAL(wp), DIMENSION(0:nx+2) ::  work   !<
748       REAL(wp), DIMENSION(nx+2)   ::  work1  !<
749       
750       COMPLEX(wp), DIMENSION(:), ALLOCATABLE ::  cwork  !<
751       
752#if defined( __ibm )
753       REAL(wp), DIMENSION(nau2) ::  aux2       !<
754       REAL(wp), DIMENSION(nau2) ::  aux4       !<
755#elif defined( __nec )
756       REAL(wp), DIMENSION(6*(nx+1)) ::  work2  !<
757#endif
758
759       IF ( direction == 'forward' )  THEN
760          forward_fft = .TRUE.
761       ELSE
762          forward_fft = .FALSE.
763       ENDIF
764
765       IF ( fft_method == 'singleton-algorithm' )  THEN
766
767!
768!--       Performing the fft with singleton's software works on every system,
769!--       since it is part of the model
770          ALLOCATE( cwork(0:nx) )
771     
772          IF ( forward_fft )   then
773
774             DO  i = 0, nx
775                cwork(i) = CMPLX( ar(i), KIND=wp )
776             ENDDO
777             ishape = SHAPE( cwork )
778             CALL FFTN( cwork, ishape )
779             DO  i = 0, (nx+1)/2
780                ar(i) = REAL( cwork(i), KIND=wp )
781             ENDDO
782             DO  i = 1, (nx+1)/2 - 1
783                ar(nx+1-i) = -AIMAG( cwork(i) )
784             ENDDO
785
786          ELSE
787
788             cwork(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
789             DO  i = 1, (nx+1)/2 - 1
790                cwork(i)      = CMPLX( ar(i), -ar(nx+1-i), KIND=wp )
791                cwork(nx+1-i) = CMPLX( ar(i),  ar(nx+1-i), KIND=wp )
792             ENDDO
793             cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp, KIND=wp )
794
795             ishape = SHAPE( cwork )
796             CALL FFTN( cwork, ishape, inv = .TRUE. )
797
798             DO  i = 0, nx
799                ar(i) = REAL( cwork(i), KIND=wp )
800             ENDDO
801
802          ENDIF
803
804          DEALLOCATE( cwork )
805
806       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
807
808!
809!--       Performing the fft with Temperton's software works on every system,
810!--       since it is part of the model
811          IF ( forward_fft )  THEN
812
813             work(0:nx) = ar
814             CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
815
816             DO  i = 0, (nx+1)/2
817                ar(i) = work(2*i)
818             ENDDO
819             DO  i = 1, (nx+1)/2 - 1
820                ar(nx+1-i) = work(2*i+1)
821             ENDDO
822
823          ELSE
824
825             DO  i = 0, (nx+1)/2
826                work(2*i) = ar(i)
827             ENDDO
828             DO  i = 1, (nx+1)/2 - 1
829                work(2*i+1) = ar(nx+1-i)
830             ENDDO
831             work(1)    = 0.0_wp
832             work(nx+2) = 0.0_wp
833
834             CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
835             ar = work(0:nx)
836
837          ENDIF
838
839       ELSEIF ( fft_method == 'fftw' )  THEN
840
841#if defined( __fftw )
842          IF ( forward_fft )  THEN
843
844             x_in(0:nx) = ar(0:nx)
845             CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
846
847             DO  i = 0, (nx+1)/2
848                ar(i) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
849             ENDDO
850             DO  i = 1, (nx+1)/2 - 1
851                ar(nx+1-i) = AIMAG( x_out(i) ) / ( nx+1 )
852             ENDDO
853
854         ELSE
855
856             x_out(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
857             DO  i = 1, (nx+1)/2 - 1
858                x_out(i) = CMPLX( ar(i), ar(nx+1-i), KIND=wp )
859             ENDDO
860             x_out((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp, KIND=wp )
861
862             CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
863             ar(0:nx) = x_in(0:nx)
864
865         ENDIF
866#endif
867
868       ELSEIF ( fft_method == 'system-specific' )  THEN
869
870#if defined( __ibm )
871          IF ( forward_fft )  THEN
872
873             CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1,   &
874                         aux2, nau2 )
875
876             DO  i = 0, (nx+1)/2
877                ar(i) = work(2*i)
878             ENDDO
879             DO  i = 1, (nx+1)/2 - 1
880                ar(nx+1-i) = work(2*i+1)
881             ENDDO
882
883          ELSE
884
885             DO  i = 0, (nx+1)/2
886                work(2*i) = ar(i)
887             ENDDO
888             DO  i = 1, (nx+1)/2 - 1
889                work(2*i+1) = ar(nx+1-i)
890             ENDDO
891             work(1) = 0.0_wp
892             work(nx+2) = 0.0_wp
893
894             CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
895                         aux4, nau2 )
896
897             DO  i = 0, nx
898                ar(i) = work(i)
899             ENDDO
900
901          ENDIF
902#elif defined( __nec )
903          IF ( forward_fft )  THEN
904
905             work(0:nx) = ar(0:nx)
906
907             CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
908     
909             DO  i = 0, (nx+1)/2
910                ar(i) = work(2*i)
911             ENDDO
912             DO  i = 1, (nx+1)/2 - 1
913                ar(nx+1-i) = work(2*i+1)
914             ENDDO
915
916          ELSE
917
918             DO  i = 0, (nx+1)/2
919                work(2*i) = ar(i)
920             ENDDO
921             DO  i = 1, (nx+1)/2 - 1
922                work(2*i+1) = ar(nx+1-i)
923             ENDDO
924             work(1) = 0.0_wp
925             work(nx+2) = 0.0_wp
926
927             CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
928
929             ar(0:nx) = work(0:nx)
930
931          ENDIF
932#else
933          message_string = 'no system-specific fft-call available'
934          CALL message( 'fft_x_1d', 'PA0188', 1, 2, 0, 6, 0 )
935#endif
936       ELSE
937          message_string = 'fft method "' // TRIM( fft_method) // &
938                           '" not available'
939          CALL message( 'fft_x_1d', 'PA0189', 1, 2, 0, 6, 0 )
940
941       ENDIF
942
943    END SUBROUTINE fft_x_1d
944
945!------------------------------------------------------------------------------!
946! Description:
947! ------------
948!> Fourier-transformation along y-direction.
949!> Version for 2D-decomposition.
950!> It uses internal algorithms (Singleton or Temperton) or
951!> system-specific routines, if they are available.
952!>
953!> direction:  'forward' or 'backward'
954!> ar, ar_tr:  3D data arrays
955!>             forward:   ar: before  ar_tr: after transformation
956!>             backward:  ar_tr: before  ar: after transfosition
957!>
958!> In case of non-overlapping transposition/transformation:
959!> nxl_y_bound = nxl_y_l = nxl_y
960!> nxr_y_bound = nxr_y_l = nxr_y
961!>
962!> In case of overlapping transposition/transformation
963!> - nxl_y_bound  and  nxr_y_bound have the original values of
964!>   nxl_y, nxr_y.  ar_tr is dimensioned using these values.
965!> - nxl_y_l = nxr_y_r.  ar is dimensioned with these values, so that
966!>   transformation is carried out for a 2D-plane only.
967!------------------------------------------------------------------------------!
968 
969    SUBROUTINE fft_y( ar, direction, ar_tr, nxl_y_bound, nxr_y_bound, nxl_y_l, &
970                      nxr_y_l )
971
972
973       IMPLICIT NONE
974
975       CHARACTER (LEN=*) ::  direction  !<
976       
977       INTEGER(iwp) ::  i            !<
978       INTEGER(iwp) ::  j            !<
979       INTEGER(iwp) ::  jshape(1)    !<
980       INTEGER(iwp) ::  k            !<
981       INTEGER(iwp) ::  nxl_y_bound  !<
982       INTEGER(iwp) ::  nxl_y_l      !<
983       INTEGER(iwp) ::  nxr_y_bound  !<
984       INTEGER(iwp) ::  nxr_y_l      !<
985
986       LOGICAL ::  forward_fft  !<
987
988       REAL(wp), DIMENSION(0:ny+2) ::  work   !<
989       REAL(wp), DIMENSION(ny+2)   ::  work1  !<
990       
991       COMPLEX(wp), DIMENSION(:), ALLOCATABLE ::  cwork  !<
992       
993#if defined( __ibm )
994       REAL(wp), DIMENSION(nau2) ::  auy2  !<
995       REAL(wp), DIMENSION(nau2) ::  auy4  !<
996#elif defined( __nec )
997       REAL(wp), DIMENSION(6*(ny+1)) ::  work2  !<
998#endif
999
1000       REAL(wp), DIMENSION(0:ny,nxl_y_l:nxr_y_l,nzb_y:nzt_y)         ::        &
1001          ar     !<
1002       REAL(wp), DIMENSION(0:ny,nxl_y_bound:nxr_y_bound,nzb_y:nzt_y) ::        &
1003          ar_tr  !<
1004
1005       IF ( direction == 'forward' )  THEN
1006          forward_fft = .TRUE.
1007       ELSE
1008          forward_fft = .FALSE.
1009       ENDIF
1010
1011       IF ( fft_method == 'singleton-algorithm' )  THEN
1012
1013!
1014!--       Performing the fft with singleton's software works on every system,
1015!--       since it is part of the model
1016          ALLOCATE( cwork(0:ny) )
1017
1018          IF ( forward_fft )   then
1019
1020             !$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
1021             !$OMP DO
1022             DO  k = nzb_y, nzt_y
1023                DO  i = nxl_y_l, nxr_y_l
1024
1025                   DO  j = 0, ny
1026                      cwork(j) = CMPLX( ar(j,i,k), KIND=wp )
1027                   ENDDO
1028
1029                   jshape = SHAPE( cwork )
1030                   CALL FFTN( cwork, jshape )
1031
1032                   DO  j = 0, (ny+1)/2
1033                      ar_tr(j,i,k) = REAL( cwork(j), KIND=wp )
1034                   ENDDO
1035                   DO  j = 1, (ny+1)/2 - 1
1036                      ar_tr(ny+1-j,i,k) = -AIMAG( cwork(j) )
1037                   ENDDO
1038
1039                ENDDO
1040             ENDDO
1041             !$OMP END PARALLEL
1042
1043          ELSE
1044
1045             !$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
1046             !$OMP DO
1047             DO  k = nzb_y, nzt_y
1048                DO  i = nxl_y_l, nxr_y_l
1049
1050                   cwork(0) = CMPLX( ar_tr(0,i,k), 0.0_wp, KIND=wp )
1051                   DO  j = 1, (ny+1)/2 - 1
1052                      cwork(j)      = CMPLX( ar_tr(j,i,k), -ar_tr(ny+1-j,i,k), &
1053                                             KIND=wp )
1054                      cwork(ny+1-j) = CMPLX( ar_tr(j,i,k),  ar_tr(ny+1-j,i,k), &
1055                                             KIND=wp )
1056                   ENDDO
1057                   cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp,       &
1058                                            KIND=wp )
1059
1060                   jshape = SHAPE( cwork )
1061                   CALL FFTN( cwork, jshape, inv = .TRUE. )
1062
1063                   DO  j = 0, ny
1064                      ar(j,i,k) = REAL( cwork(j), KIND=wp )
1065                   ENDDO
1066
1067                ENDDO
1068             ENDDO
1069             !$OMP END PARALLEL
1070
1071          ENDIF
1072
1073          DEALLOCATE( cwork )
1074
1075       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
1076
1077!
1078!--       Performing the fft with Temperton's software works on every system,
1079!--       since it is part of the model
1080          IF ( forward_fft )  THEN
1081
1082             !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
1083             !$OMP DO
1084             DO  k = nzb_y, nzt_y
1085                DO  i = nxl_y_l, nxr_y_l
1086
1087                   work(0:ny) = ar(0:ny,i,k)
1088                   CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
1089
1090                   DO  j = 0, (ny+1)/2
1091                      ar_tr(j,i,k) = work(2*j)
1092                   ENDDO
1093                   DO  j = 1, (ny+1)/2 - 1
1094                      ar_tr(ny+1-j,i,k) = work(2*j+1)
1095                   ENDDO
1096
1097                ENDDO
1098             ENDDO
1099             !$OMP END PARALLEL
1100
1101          ELSE
1102
1103             !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
1104             !$OMP DO
1105             DO  k = nzb_y, nzt_y
1106                DO  i = nxl_y_l, nxr_y_l
1107
1108                   DO  j = 0, (ny+1)/2
1109                      work(2*j) = ar_tr(j,i,k)
1110                   ENDDO
1111                   DO  j = 1, (ny+1)/2 - 1
1112                      work(2*j+1) = ar_tr(ny+1-j,i,k)
1113                   ENDDO
1114                   work(1)    = 0.0_wp
1115                   work(ny+2) = 0.0_wp
1116
1117                   CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
1118                   ar(0:ny,i,k) = work(0:ny)
1119
1120                ENDDO
1121             ENDDO
1122             !$OMP END PARALLEL
1123
1124          ENDIF
1125
1126       ELSEIF ( fft_method == 'fftw' )  THEN
1127
1128#if defined( __fftw )
1129          IF ( forward_fft )  THEN
1130
1131             !$OMP PARALLEL PRIVATE ( work, i, j, k )
1132             !$OMP DO
1133             DO  k = nzb_y, nzt_y
1134                DO  i = nxl_y_l, nxr_y_l
1135
1136                   y_in(0:ny) = ar(0:ny,i,k)
1137                   CALL FFTW_EXECUTE_DFT_R2C( plan_yf, y_in, y_out )
1138
1139                   DO  j = 0, (ny+1)/2
1140                      ar_tr(j,i,k) = REAL( y_out(j), KIND=wp ) / (ny+1)
1141                   ENDDO
1142                   DO  j = 1, (ny+1)/2 - 1
1143                      ar_tr(ny+1-j,i,k) = AIMAG( y_out(j) ) / (ny+1)
1144                   ENDDO
1145
1146                ENDDO
1147             ENDDO
1148             !$OMP END PARALLEL
1149
1150          ELSE
1151
1152             !$OMP PARALLEL PRIVATE ( work, i, j, k )
1153             !$OMP DO
1154             DO  k = nzb_y, nzt_y
1155                DO  i = nxl_y_l, nxr_y_l
1156
1157                   y_out(0) = CMPLX( ar_tr(0,i,k), 0.0_wp, KIND=wp )
1158                   DO  j = 1, (ny+1)/2 - 1
1159                      y_out(j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k),       &
1160                                        KIND=wp )
1161                   ENDDO
1162                   y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp,       &
1163                                            KIND=wp )
1164
1165                   CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
1166                   ar(0:ny,i,k) = y_in(0:ny)
1167
1168                ENDDO
1169             ENDDO
1170             !$OMP END PARALLEL
1171
1172          ENDIF
1173#endif
1174
1175       ELSEIF ( fft_method == 'system-specific' )  THEN
1176
1177#if defined( __ibm )
1178          IF ( forward_fft)  THEN
1179
1180             !$OMP PARALLEL PRIVATE ( work, i, j, k )
1181             !$OMP DO
1182             DO  k = nzb_y, nzt_y
1183                DO  i = nxl_y_l, nxr_y_l
1184
1185                   CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1,   & 
1186                               nau1, auy2, nau2 )
1187
1188                   DO  j = 0, (ny+1)/2
1189                      ar_tr(j,i,k) = work(2*j)
1190                   ENDDO
1191                   DO  j = 1, (ny+1)/2 - 1
1192                      ar_tr(ny+1-j,i,k) = work(2*j+1)
1193                   ENDDO
1194
1195                ENDDO
1196             ENDDO
1197             !$OMP END PARALLEL
1198
1199          ELSE
1200
1201             !$OMP PARALLEL PRIVATE ( work, i, j, k )
1202             !$OMP DO
1203             DO  k = nzb_y, nzt_y
1204                DO  i = nxl_y_l, nxr_y_l
1205
1206                   DO  j = 0, (ny+1)/2
1207                      work(2*j) = ar_tr(j,i,k)
1208                   ENDDO
1209                   DO  j = 1, (ny+1)/2 - 1
1210                      work(2*j+1) = ar_tr(ny+1-j,i,k)
1211                   ENDDO
1212                   work(1)    = 0.0_wp
1213                   work(ny+2) = 0.0_wp
1214
1215                   CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny,      &
1216                               auy3, nau1, auy4, nau2 )
1217
1218                   DO  j = 0, ny
1219                      ar(j,i,k) = work(j)
1220                   ENDDO
1221
1222                ENDDO
1223             ENDDO
1224             !$OMP END PARALLEL
1225
1226          ENDIF
1227#elif defined( __nec )
1228          IF ( forward_fft )  THEN
1229
1230             !$OMP PARALLEL PRIVATE ( work, i, j, k )
1231             !$OMP DO
1232             DO  k = nzb_y, nzt_y
1233                DO  i = nxl_y_l, nxr_y_l
1234
1235                   work(0:ny) = ar(0:ny,i,k)
1236
1237                   CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
1238
1239                   DO  j = 0, (ny+1)/2
1240                      ar_tr(j,i,k) = work(2*j)
1241                   ENDDO
1242                   DO  j = 1, (ny+1)/2 - 1
1243                      ar_tr(ny+1-j,i,k) = work(2*j+1)
1244                   ENDDO
1245
1246                ENDDO
1247             ENDDO
1248             !$END OMP PARALLEL
1249
1250          ELSE
1251
1252             !$OMP PARALLEL PRIVATE ( work, i, j, k )
1253             !$OMP DO
1254             DO  k = nzb_y, nzt_y
1255                DO  i = nxl_y_l, nxr_y_l
1256
1257                   DO  j = 0, (ny+1)/2
1258                      work(2*j) = ar_tr(j,i,k)
1259                   ENDDO
1260                   DO  j = 1, (ny+1)/2 - 1
1261                      work(2*j+1) = ar_tr(ny+1-j,i,k)
1262                   ENDDO
1263                   work(1) = 0.0_wp
1264                   work(ny+2) = 0.0_wp
1265
1266                   CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1267
1268                   ar(0:ny,i,k) = work(0:ny)
1269
1270                ENDDO
1271             ENDDO
1272             !$OMP END PARALLEL
1273
1274          ENDIF
1275#else
1276          message_string = 'no system-specific fft-call available'
1277          CALL message( 'fft_y', 'PA0188', 1, 2, 0, 6, 0 ) 
1278#endif
1279
1280       ELSE
1281
1282          message_string = 'fft method "' // TRIM( fft_method) // &
1283                           '" not available'
1284          CALL message( 'fft_y', 'PA0189', 1, 2, 0, 6, 0 )
1285
1286       ENDIF
1287
1288    END SUBROUTINE fft_y
1289
1290!------------------------------------------------------------------------------!
1291! Description:
1292! ------------
1293!> Fourier-transformation along y-direction.
1294!> Version for 1D-decomposition.
1295!> It uses internal algorithms (Singleton or Temperton) or
1296!> system-specific routines, if they are available.
1297!------------------------------------------------------------------------------!
1298 
1299    SUBROUTINE fft_y_1d( ar, direction )
1300
1301
1302       IMPLICIT NONE
1303
1304       CHARACTER (LEN=*) ::  direction
1305       
1306       INTEGER(iwp) ::  j          !<
1307       INTEGER(iwp) ::  jshape(1)  !<
1308
1309       LOGICAL ::  forward_fft  !<
1310
1311       REAL(wp), DIMENSION(0:ny)    ::  ar     !<
1312       REAL(wp), DIMENSION(0:ny+2)  ::  work   !<
1313       REAL(wp), DIMENSION(ny+2)    ::  work1  !<
1314       
1315       COMPLEX(wp), DIMENSION(:), ALLOCATABLE ::  cwork  !<
1316       
1317#if defined( __ibm )
1318       REAL(wp), DIMENSION(nau2) ::  auy2  !<
1319       REAL(wp), DIMENSION(nau2) ::  auy4  !<
1320#elif defined( __nec )
1321       REAL(wp), DIMENSION(6*(ny+1)) ::  work2  !<
1322#endif
1323
1324       IF ( direction == 'forward' )  THEN
1325          forward_fft = .TRUE.
1326       ELSE
1327          forward_fft = .FALSE.
1328       ENDIF
1329
1330       IF ( fft_method == 'singleton-algorithm' )  THEN
1331
1332!
1333!--       Performing the fft with singleton's software works on every system,
1334!--       since it is part of the model
1335          ALLOCATE( cwork(0:ny) )
1336
1337          IF ( forward_fft )  THEN
1338
1339             DO  j = 0, ny
1340                cwork(j) = CMPLX( ar(j), KIND=wp )
1341             ENDDO
1342
1343             jshape = SHAPE( cwork )
1344             CALL FFTN( cwork, jshape )
1345
1346             DO  j = 0, (ny+1)/2
1347                ar(j) = REAL( cwork(j), KIND=wp )
1348             ENDDO
1349             DO  j = 1, (ny+1)/2 - 1
1350                ar(ny+1-j) = -AIMAG( cwork(j) )
1351             ENDDO
1352
1353          ELSE
1354
1355             cwork(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
1356             DO  j = 1, (ny+1)/2 - 1
1357                cwork(j)      = CMPLX( ar(j), -ar(ny+1-j), KIND=wp )
1358                cwork(ny+1-j) = CMPLX( ar(j),  ar(ny+1-j), KIND=wp )
1359             ENDDO
1360             cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp, KIND=wp )
1361
1362             jshape = SHAPE( cwork )
1363             CALL FFTN( cwork, jshape, inv = .TRUE. )
1364
1365             DO  j = 0, ny
1366                ar(j) = REAL( cwork(j), KIND=wp )
1367             ENDDO
1368
1369          ENDIF
1370
1371          DEALLOCATE( cwork )
1372
1373       ELSEIF ( fft_method == 'temperton-algorithm' )  THEN
1374
1375!
1376!--       Performing the fft with Temperton's software works on every system,
1377!--       since it is part of the model
1378          IF ( forward_fft )  THEN
1379
1380             work(0:ny) = ar
1381             CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
1382
1383             DO  j = 0, (ny+1)/2
1384                ar(j) = work(2*j)
1385             ENDDO
1386             DO  j = 1, (ny+1)/2 - 1
1387                ar(ny+1-j) = work(2*j+1)
1388             ENDDO
1389
1390          ELSE
1391
1392             DO  j = 0, (ny+1)/2
1393                work(2*j) = ar(j)
1394             ENDDO
1395             DO  j = 1, (ny+1)/2 - 1
1396                work(2*j+1) = ar(ny+1-j)
1397             ENDDO
1398             work(1)    = 0.0_wp
1399             work(ny+2) = 0.0_wp
1400
1401             CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
1402             ar = work(0:ny)
1403
1404          ENDIF
1405
1406       ELSEIF ( fft_method == 'fftw' )  THEN
1407
1408#if defined( __fftw )
1409          IF ( forward_fft )  THEN
1410
1411             y_in(0:ny) = ar(0:ny)
1412             CALL FFTW_EXECUTE_DFT_R2C( plan_yf, y_in, y_out )
1413
1414             DO  j = 0, (ny+1)/2
1415                ar(j) = REAL( y_out(j), KIND=wp ) / (ny+1)
1416             ENDDO
1417             DO  j = 1, (ny+1)/2 - 1
1418                ar(ny+1-j) = AIMAG( y_out(j) ) / (ny+1)
1419             ENDDO
1420
1421          ELSE
1422
1423             y_out(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
1424             DO  j = 1, (ny+1)/2 - 1
1425                y_out(j) = CMPLX( ar(j), ar(ny+1-j), KIND=wp )
1426             ENDDO
1427             y_out((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp, KIND=wp )
1428
1429             CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
1430             ar(0:ny) = y_in(0:ny)
1431
1432          ENDIF
1433#endif
1434
1435       ELSEIF ( fft_method == 'system-specific' )  THEN
1436
1437#if defined( __ibm )
1438          IF ( forward_fft )  THEN
1439
1440             CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1,   &
1441                         auy2, nau2 )
1442
1443             DO  j = 0, (ny+1)/2
1444                ar(j) = work(2*j)
1445             ENDDO
1446             DO  j = 1, (ny+1)/2 - 1
1447                ar(ny+1-j) = work(2*j+1)
1448             ENDDO
1449
1450          ELSE
1451
1452             DO  j = 0, (ny+1)/2
1453                work(2*j) = ar(j)
1454             ENDDO
1455             DO  j = 1, (ny+1)/2 - 1
1456                work(2*j+1) = ar(ny+1-j)
1457             ENDDO
1458             work(1)    = 0.0_wp
1459             work(ny+2) = 0.0_wp
1460
1461             CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3,      &
1462                         nau1, auy4, nau2 )
1463
1464             DO  j = 0, ny
1465                ar(j) = work(j)
1466             ENDDO
1467
1468          ENDIF
1469#elif defined( __nec )
1470          IF ( forward_fft )  THEN
1471
1472             work(0:ny) = ar(0:ny)
1473
1474             CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
1475
1476             DO  j = 0, (ny+1)/2
1477                ar(j) = work(2*j)
1478             ENDDO
1479             DO  j = 1, (ny+1)/2 - 1
1480                ar(ny+1-j) = work(2*j+1)
1481             ENDDO
1482
1483          ELSE
1484
1485             DO  j = 0, (ny+1)/2
1486                work(2*j) = ar(j)
1487             ENDDO
1488             DO  j = 1, (ny+1)/2 - 1
1489                work(2*j+1) = ar(ny+1-j)
1490             ENDDO
1491             work(1) = 0.0_wp
1492             work(ny+2) = 0.0_wp
1493
1494             CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
1495
1496             ar(0:ny) = work(0:ny)
1497
1498          ENDIF
1499#else
1500          message_string = 'no system-specific fft-call available'
1501          CALL message( 'fft_y_1d', 'PA0188', 1, 2, 0, 6, 0 ) 
1502
1503#endif
1504
1505       ELSE
1506
1507          message_string = 'fft method "' // TRIM( fft_method) // &
1508                           '" not available'
1509          CALL message( 'fft_y_1d', 'PA0189', 1, 2, 0, 6, 0 )
1510
1511       ENDIF
1512
1513    END SUBROUTINE fft_y_1d
1514
1515!------------------------------------------------------------------------------!
1516! Description:
1517! ------------
1518!> Fourier-transformation along x-direction.
1519!> Version for 1d domain decomposition
1520!> using multiple 1D FFT from Math Keisan on NEC or Temperton-algorithm
1521!> (no singleton-algorithm on NEC because it does not vectorize)
1522!------------------------------------------------------------------------------!
1523 
1524    SUBROUTINE fft_x_m( ar, direction )
1525
1526
1527       IMPLICIT NONE
1528
1529       CHARACTER (LEN=*) ::  direction  !<
1530       
1531       INTEGER(iwp) ::  i     !<
1532       INTEGER(iwp) ::  k     !<
1533       INTEGER(iwp) ::  siza  !<
1534
1535       REAL(wp), DIMENSION(0:nx,nz)       ::  ar     !<
1536       REAL(wp), DIMENSION(0:nx+3,nz+1)   ::  ai     !<
1537       REAL(wp), DIMENSION(6*(nx+4),nz+1) ::  work1  !<
1538       
1539#if defined( __nec )
1540       INTEGER(iwp) ::  sizw  !<
1541       
1542       COMPLEX(wp), DIMENSION((nx+4)/2+1,nz+1) ::  work  !<
1543#endif
1544
1545       IF ( fft_method == 'temperton-algorithm' )  THEN
1546
1547          siza = SIZE( ai, 1 )
1548
1549          IF ( direction == 'forward')  THEN
1550
1551             ai(0:nx,1:nz) = ar(0:nx,1:nz)
1552             ai(nx+1:,:)   = 0.0_wp
1553
1554             CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
1555
1556             DO  k = 1, nz
1557                DO  i = 0, (nx+1)/2
1558                   ar(i,k) = ai(2*i,k)
1559                ENDDO
1560                DO  i = 1, (nx+1)/2 - 1
1561                   ar(nx+1-i,k) = ai(2*i+1,k)
1562                ENDDO
1563             ENDDO
1564
1565          ELSE
1566
1567             DO  k = 1, nz
1568                DO  i = 0, (nx+1)/2
1569                   ai(2*i,k) = ar(i,k)
1570                ENDDO
1571                DO  i = 1, (nx+1)/2 - 1
1572                   ai(2*i+1,k) = ar(nx+1-i,k)
1573                ENDDO
1574                ai(1,k) = 0.0_wp
1575                ai(nx+2,k) = 0.0_wp
1576             ENDDO
1577
1578             CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
1579
1580             ar(0:nx,1:nz) = ai(0:nx,1:nz)
1581
1582          ENDIF
1583
1584       ELSEIF ( fft_method == 'system-specific' )  THEN
1585
1586#if defined( __nec )
1587          siza = SIZE( ai, 1 )
1588          sizw = SIZE( work, 1 )
1589
1590          IF ( direction == 'forward')  THEN
1591
1592!
1593!--          Tables are initialized once more. This call should not be
1594!--          necessary, but otherwise program aborts in asymmetric case
1595             CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4,       &
1596                          trig_xf, work1, 0 )
1597
1598             ai(0:nx,1:nz) = ar(0:nx,1:nz)
1599             IF ( nz1 > nz )  THEN
1600                ai(:,nz1) = 0.0_wp
1601             ENDIF
1602
1603             CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw,         &
1604                          trig_xf, work1, 0 )
1605
1606             DO  k = 1, nz
1607                DO  i = 0, (nx+1)/2
1608                   ar(i,k) = REAL( work(i+1,k), KIND=wp )
1609                ENDDO
1610                DO  i = 1, (nx+1)/2 - 1
1611                   ar(nx+1-i,k) = AIMAG( work(i+1,k) )
1612                ENDDO
1613             ENDDO
1614
1615          ELSE
1616
1617!
1618!--          Tables are initialized once more. This call should not be
1619!--          necessary, but otherwise program aborts in asymmetric case
1620             CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4,       &
1621                          trig_xb, work1, 0 )
1622
1623             IF ( nz1 > nz )  THEN
1624                work(:,nz1) = 0.0_wp
1625             ENDIF
1626             DO  k = 1, nz
1627                work(1,k) = CMPLX( ar(0,k), 0.0_wp, KIND=wp )
1628                DO  i = 1, (nx+1)/2 - 1
1629                   work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k), KIND=wp )
1630                ENDDO
1631                work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0_wp, KIND=wp )
1632             ENDDO
1633
1634             CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, &
1635                          trig_xb, work1, 0 )
1636
1637             ar(0:nx,1:nz) = ai(0:nx,1:nz)
1638
1639          ENDIF
1640
1641#else
1642          message_string = 'no system-specific fft-call available'
1643          CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 ) 
1644#endif
1645
1646       ELSE
1647
1648          message_string = 'fft method "' // TRIM( fft_method) // &
1649                           '" not available'
1650          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
1651
1652       ENDIF
1653
1654    END SUBROUTINE fft_x_m
1655
1656!------------------------------------------------------------------------------!
1657! Description:
1658! ------------
1659!> Fourier-transformation along y-direction.
1660!> Version for 1d domain decomposition
1661!> using multiple 1D FFT from Math Keisan on NEC or Temperton-algorithm
1662!> (no singleton-algorithm on NEC because it does not vectorize)
1663!------------------------------------------------------------------------------!
1664 
1665    SUBROUTINE fft_y_m( ar, ny1, direction )
1666
1667
1668       IMPLICIT NONE
1669
1670       CHARACTER (LEN=*) ::  direction  !<
1671       
1672       INTEGER(iwp) ::  j     !<
1673       INTEGER(iwp) ::  k     !<
1674       INTEGER(iwp) ::  ny1   !<
1675       INTEGER(iwp) ::  siza  !<
1676
1677       REAL(wp), DIMENSION(0:ny1,nz)      ::  ar     !<
1678       REAL(wp), DIMENSION(0:ny+3,nz+1)   ::  ai     !<
1679       REAL(wp), DIMENSION(6*(ny+4),nz+1) ::  work1  !<
1680       
1681#if defined( __nec )
1682       INTEGER(iwp) ::  sizw  !<
1683       
1684       COMPLEX(wp), DIMENSION((ny+4)/2+1,nz+1) ::  work !<
1685#endif
1686
1687       IF ( fft_method == 'temperton-algorithm' )  THEN
1688
1689          siza = SIZE( ai, 1 )
1690
1691          IF ( direction == 'forward')  THEN
1692
1693             ai(0:ny,1:nz) = ar(0:ny,1:nz)
1694             ai(ny+1:,:)   = 0.0_wp
1695
1696             CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
1697
1698             DO  k = 1, nz
1699                DO  j = 0, (ny+1)/2
1700                   ar(j,k) = ai(2*j,k)
1701                ENDDO
1702                DO  j = 1, (ny+1)/2 - 1
1703                   ar(ny+1-j,k) = ai(2*j+1,k)
1704                ENDDO
1705             ENDDO
1706
1707          ELSE
1708
1709             DO  k = 1, nz
1710                DO  j = 0, (ny+1)/2
1711                   ai(2*j,k) = ar(j,k)
1712                ENDDO
1713                DO  j = 1, (ny+1)/2 - 1
1714                   ai(2*j+1,k) = ar(ny+1-j,k)
1715                ENDDO
1716                ai(1,k) = 0.0_wp
1717                ai(ny+2,k) = 0.0_wp
1718             ENDDO
1719
1720             CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
1721
1722             ar(0:ny,1:nz) = ai(0:ny,1:nz)
1723
1724          ENDIF
1725
1726       ELSEIF ( fft_method == 'system-specific' )  THEN
1727
1728#if defined( __nec )
1729          siza = SIZE( ai, 1 )
1730          sizw = SIZE( work, 1 )
1731
1732          IF ( direction == 'forward')  THEN
1733
1734!
1735!--          Tables are initialized once more. This call should not be
1736!--          necessary, but otherwise program aborts in asymmetric case
1737             CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
1738                          trig_yf, work1, 0 )
1739
1740             ai(0:ny,1:nz) = ar(0:ny,1:nz)
1741             IF ( nz1 > nz )  THEN
1742                ai(:,nz1) = 0.0_wp
1743             ENDIF
1744
1745             CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, &
1746                          trig_yf, work1, 0 )
1747
1748             DO  k = 1, nz
1749                DO  j = 0, (ny+1)/2
1750                   ar(j,k) = REAL( work(j+1,k), KIND=wp )
1751                ENDDO
1752                DO  j = 1, (ny+1)/2 - 1
1753                   ar(ny+1-j,k) = AIMAG( work(j+1,k) )
1754                ENDDO
1755             ENDDO
1756
1757          ELSE
1758
1759!
1760!--          Tables are initialized once more. This call should not be
1761!--          necessary, but otherwise program aborts in asymmetric case
1762             CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
1763                          trig_yb, work1, 0 )
1764
1765             IF ( nz1 > nz )  THEN
1766                work(:,nz1) = 0.0_wp
1767             ENDIF
1768             DO  k = 1, nz
1769                work(1,k) = CMPLX( ar(0,k), 0.0_wp, KIND=wp )
1770                DO  j = 1, (ny+1)/2 - 1
1771                   work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k), KIND=wp )
1772                ENDDO
1773                work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0_wp, KIND=wp )
1774             ENDDO
1775
1776             CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, &
1777                          trig_yb, work1, 0 )
1778
1779             ar(0:ny,1:nz) = ai(0:ny,1:nz)
1780
1781          ENDIF
1782
1783#else
1784          message_string = 'no system-specific fft-call available'
1785          CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 ) 
1786#endif
1787
1788       ELSE
1789         
1790          message_string = 'fft method "' // TRIM( fft_method) // &
1791                           '" not available'
1792          CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
1793
1794       ENDIF
1795
1796    END SUBROUTINE fft_y_m
1797
1798
1799 END MODULE fft_xy
Note: See TracBrowser for help on using the repository browser.