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

Last change on this file since 2287 was 2274, checked in by Giersch, 7 years ago

Changed error messages

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