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

Last change on this file since 3560 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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