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

Last change on this file since 3045 was 3045, checked in by Giersch, 3 years ago

Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled

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