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

Last change on this file since 2258 was 2119, checked in by raasch, 8 years ago

last commit documented

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