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

Last change on this file since 2716 was 2716, checked in by kanani, 6 years ago

Correction of "Former revisions" section

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