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

Last change on this file since 2416 was 2300, checked in by raasch, 7 years ago

NEC related code partly removed, host variable partly removed, host specific code completely removed, default values for host, loop_optimization and termination time_needed changed

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