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

Last change on this file since 3600 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
RevLine 
[1850]1!> @file fft_xy_mod.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]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.
[1036]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!
[2718]17! Copyright 1997-2018 Leibniz Universitaet Hannover
[1322]18!------------------------------------------------------------------------------!
[1036]19!
[254]20! Current revisions:
[1]21! -----------------
[1683]22!
[2119]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: fft_xy_mod.f90 3241 2018-09-12 15:02:00Z banzhafs $
[3241]27! preprocessor switches for variables that are required on NEC only
28!
29! 3045 2018-05-28 07:55:41Z Giersch
[3045]30! Error messages revised
31!
32! 2718 2018-01-02 08:49:38Z maronga
[2716]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
[2300]39! NEC related code partly removed, host replaced by loop_optimization
40!
41! 2274 2017-06-09 13:27:48Z Giersch
[2274]42! Changed error messages
43!
44! 2119 2017-01-17 16:51:50Z raasch
[1321]45!
[2119]46! 2118 2017-01-17 16:38:49Z raasch
47! OpenACC directives and CUDA-fft related code removed
48!
[2001]49! 2000 2016-08-20 18:09:15Z knoop
50! Forced header and separation lines into 80 columns
51!
[1851]52! 1850 2016-04-08 13:29:27Z maronga
53! Module renamed
54!
[1816]55! 1815 2016-04-06 13:49:59Z raasch
56! cpp-directives for ibmy removed
57!
[1750]58! 1749 2016-02-09 12:19:56Z raasch
59! small OpenACC bugfix
60!
[1683]61! 1682 2015-10-07 23:56:08Z knoop
62! Code annotations made doxygen readable
63!
[1601]64! 1600 2015-06-11 15:50:12Z raasch
65! bugfix: openMP threadprivate statement moved after variable declaration
66!
[1483]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!
[1403]70! 1402 2014-05-09 14:25:13Z raasch
71! fortran bugfix for r1392
72!
[1399]73! 1398 2014-05-07 11:15:00Z heinze
74! bugfix: typo removed for KIND in CMPLX function
75!
[1393]76! 1392 2014-05-06 09:10:05Z raasch
77! bugfix: KIND attribute added to CMPLX functions
78!
[1375]79! 1374 2014-04-25 12:55:07Z raasch
80! bugfixes: missing variables added to ONLY list, dpk renamed dp
81!
[1373]82! 1372 2014-04-24 06:29:32Z raasch
83! openMP-bugfix for fftw: some arrays defined as threadprivate
84!
[1354]85! 1353 2014-04-08 15:21:23Z heinze
86! REAL constants provided with KIND-attribute
87!
[1343]88! 1342 2014-03-26 17:04:47Z kanani
89! REAL constants defined as wp-kind
90!
[1323]91! 1322 2014-03-20 16:38:49Z raasch
92! REAL functions provided with KIND-attribute
93!
[1321]94! 1320 2014-03-20 08:40:49Z raasch
[1320]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
[1]102!
[1305]103! 1304 2014-03-12 10:29:42Z raasch
104! openmp bugfix: work1 used in Temperton algorithm must be private
105!
[1258]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!
[1220]110! 1219 2013-08-30 09:33:18Z heinze
111! bugfix: use own branch for fftw
112!
[1217]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!
[1211]118! 1210 2013-08-14 10:58:20Z raasch
119! fftw added
120!
[1167]121! 1166 2013-05-24 13:55:44Z raasch
122! C_DOUBLE/COMPLEX reset to dpk
123!
[1154]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!
[1112]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!
[1107]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!
[1093]137! 1092 2013-02-02 11:24:22Z raasch
138! variable sizw declared for NEC case only
139!
[1037]140! 1036 2012-10-22 13:43:42Z raasch
141! code put under GPL (PALM 3.9)
142!
[1]143! Revision 1.1  2002/06/11 13:00:49  raasch
144! Initial revision
145!
146!
147! Description:
148! ------------
[1682]149!> Fast Fourier transformation along x and y for 1d domain decomposition along x.
150!> Original version: Klaus Ketelsen (May 2002)
[1]151!------------------------------------------------------------------------------!
[1682]152 MODULE fft_xy
153 
[1]154
[1320]155    USE control_parameters,                                                    &
156        ONLY:  fft_method, message_string
157       
158    USE indices,                                                               &
159        ONLY:  nx, ny, nz
160       
[2118]161#if defined( __fftw )
[1210]162    USE, INTRINSIC ::  ISO_C_BINDING
[1153]163#endif
[1320]164
165    USE kinds
166   
167    USE singleton,                                                             &
168        ONLY: fftn
169   
[1]170    USE temperton_fft
[1320]171   
172    USE transpose_indices,                                                     &
[1374]173        ONLY:  nxl_y, nxr_y, nyn_x, nys_x, nzb_x, nzb_y, nzt_x, nzt_y
[1]174
175    IMPLICIT NONE
176
177    PRIVATE
[1106]178    PUBLIC fft_x, fft_x_1d, fft_y, fft_y_1d, fft_init, fft_x_m, fft_y_m
[1]179
[1682]180    INTEGER(iwp), DIMENSION(:), ALLOCATABLE, SAVE ::  ifax_x  !<
181    INTEGER(iwp), DIMENSION(:), ALLOCATABLE, SAVE ::  ifax_y  !<
[1]182
[1682]183    LOGICAL, SAVE ::  init_fft = .FALSE.  !<
[1]184
[1682]185    REAL(wp), SAVE ::  dnx      !<
186    REAL(wp), SAVE ::  dny      !<
187    REAL(wp), SAVE ::  sqr_dnx  !<
188    REAL(wp), SAVE ::  sqr_dny  !<
[1320]189   
[1682]190    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trigs_x  !<
191    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trigs_y  !<
[1]192
193#if defined( __ibm )
[1682]194    INTEGER(iwp), PARAMETER ::  nau1 = 20000  !<
195    INTEGER(iwp), PARAMETER ::  nau2 = 22000  !<
[1]196!
197!-- The following working arrays contain tables and have to be "save" and
198!-- shared in OpenMP sense
[1682]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  !<
[1320]203   
[1]204#elif defined( __nec )
[1682]205    INTEGER(iwp), SAVE ::  nz1  !<
[1320]206   
[1682]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  !<
[1320]211   
[1219]212#endif
213
214#if defined( __fftw )
[1210]215    INCLUDE  'fftw3.f03'
[1682]216    INTEGER(KIND=C_INT) ::  nx_c  !<
217    INTEGER(KIND=C_INT) ::  ny_c  !<
[1320]218   
[1682]219    COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE ::  x_out  !<
[1320]220    COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE ::         &
[1682]221       y_out  !<
[1320]222   
223    REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE ::                    &
[1682]224       x_in   !<
[1320]225    REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE ::                    &
[1682]226       y_in   !<
[1600]227    !$OMP THREADPRIVATE( x_out, y_out, x_in, y_in )
[1320]228   
229   
[1210]230    TYPE(C_PTR), SAVE ::  plan_xf, plan_xi, plan_yf, plan_yi
[1]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
[1106]243    INTERFACE fft_x_1d
244       MODULE PROCEDURE fft_x_1d
245    END INTERFACE fft_x_1d
246
[1]247    INTERFACE fft_y
248       MODULE PROCEDURE fft_y
249    END INTERFACE fft_y
250
[1106]251    INTERFACE fft_y_1d
252       MODULE PROCEDURE fft_y_1d
253    END INTERFACE fft_y_1d
254
[1]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
[1682]266!------------------------------------------------------------------------------!
267! Description:
268! ------------
269!> @todo Missing subroutine description.
270!------------------------------------------------------------------------------!
[1]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 )
[1682]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   !<
[1]285#elif defined( __nec )
[1682]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   !<
[1]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
[1342]302          dnx = 1.0_wp / ( nx + 1.0_wp )
303          dny = 1.0_wp / ( ny + 1.0_wp )
[1106]304          sqr_dnx = SQRT( dnx )
305          sqr_dny = SQRT( dny )
[1815]306#if defined( __ibm )
[1]307!
308!--       Initialize tables for fft along x
[1106]309          CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1,  1, sqr_dnx, aux1, nau1, &
[1]310                      aux2, nau2 )
[1106]311          CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
[1]312                      aux4, nau2 )
313!
314!--       Initialize tables for fft along y
[1106]315          CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1,  1, sqr_dny, auy1, nau1, &
[1]316                      auy2, nau2 )
[1106]317          CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
[1]318                      auy4, nau2 )
319#elif defined( __nec )
[254]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 )
[1]323
[1320]324          ALLOCATE( trig_xb(2*(nx+1)), trig_xf(2*(nx+1)),                      &
[1]325                    trig_yb(2*(ny+1)), trig_yf(2*(ny+1)) )
326
[1342]327          work_x = 0.0_wp
328          work_y = 0.0_wp
[1]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))
[1106]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 )
[1320]336          CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4,      &
[1]337                       trig_xf, workx, 0 )
[1320]338          CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4,      &
[1]339                       trig_xb, workx, 0 )
340!
341!--       Initialize tables for fft along y (non-vector and vector case (M))
[1106]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 )
[1320]344          CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4,      &
[1]345                       trig_yf, worky, 0 )
[1320]346          CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4,      &
[1]347                       trig_yb, worky, 0 )
348#else
[254]349          message_string = 'no system-specific fft-call available'
350          CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
[1]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
[1210]361       ELSEIF ( fft_method == 'fftw' )  THEN
362!
363!--       FFTW
364#if defined( __fftw )
365          nx_c = nx+1
366          ny_c = ny+1
[1372]367          !$OMP PARALLEL
[1320]368          ALLOCATE( x_in(0:nx+2), y_in(0:ny+2), x_out(0:(nx+1)/2),             &
[1210]369                    y_out(0:(ny+1)/2) )
[1372]370          !$OMP END PARALLEL
[1210]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
[1]380       ELSEIF ( fft_method == 'singleton-algorithm' )  THEN
381
382          CONTINUE
383
384       ELSE
385
[254]386          message_string = 'fft method "' // TRIM( fft_method) // &
387                           '" not available'
388          CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
[1]389       ENDIF
390
391    END SUBROUTINE fft_init
392
393
[1682]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 
[1216]403    SUBROUTINE fft_x( ar, direction, ar_2d )
[1]404
405
406       IMPLICIT NONE
407
[1682]408       CHARACTER (LEN=*) ::  direction  !<
[1320]409       
[1682]410       COMPLEX(wp), DIMENSION(:), ALLOCATABLE ::  cwork  !<
[1106]411
[1682]412       INTEGER(iwp) ::  i          !<
413       INTEGER(iwp) ::  ishape(1)  !<
414       INTEGER(iwp) ::  j          !<
415       INTEGER(iwp) ::  k          !<
[1106]416
[1682]417       LOGICAL ::  forward_fft !<
[1320]418       
[1682]419       REAL(wp), DIMENSION(0:nx+2) ::  work   !<
420       REAL(wp), DIMENSION(nx+2)   ::  work1  !<
[1320]421       
[1106]422#if defined( __ibm )
[1682]423       REAL(wp), DIMENSION(nau2) ::  aux2  !<
424       REAL(wp), DIMENSION(nau2) ::  aux4  !<
[1106]425#elif defined( __nec )
[1682]426       REAL(wp), DIMENSION(6*(nx+1)) ::  work2  !<
[1106]427#endif
428
[1320]429       REAL(wp), DIMENSION(0:nx,nys_x:nyn_x), OPTIONAL   ::                    &
[1682]430          ar_2d   !<
[1320]431       REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::                    &
[1682]432          ar      !<
[1320]433
[1106]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
[1392]455                      cwork(i) = CMPLX( ar(i,j,k), KIND=wp )
[1106]456                   ENDDO
457
458                   ishape = SHAPE( cwork )
459                   CALL FFTN( cwork, ishape )
460
461                   DO  i = 0, (nx+1)/2
[1322]462                      ar(i,j,k) = REAL( cwork(i), KIND=wp )
[1106]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
[1392]479                   cwork(0) = CMPLX( ar(0,j,k), 0.0_wp, KIND=wp )
[1106]480                   DO  i = 1, (nx+1)/2 - 1
[1392]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 )
[1106]485                   ENDDO
[1392]486                   cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp, KIND=wp )
[1106]487
488                   ishape = SHAPE( cwork )
489                   CALL FFTN( cwork, ishape, inv = .TRUE. )
490
491                   DO  i = 0, nx
[1322]492                      ar(i,j,k) = REAL( cwork(i), KIND=wp )
[1106]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
[1304]510             !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
[1106]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
[1304]531             !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
[1106]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
[1342]542                   work(1)    = 0.0_wp
543                   work(nx+2) = 0.0_wp
[1106]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
[1210]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
[1216]567                   IF ( PRESENT( ar_2d ) )  THEN
[1210]568
[1216]569                      DO  i = 0, (nx+1)/2
[1322]570                         ar_2d(i,j) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
[1216]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
[1322]579                         ar(i,j,k) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
[1216]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
[1210]587                ENDDO
588             ENDDO
589             !$OMP END PARALLEL
590
[1216]591          ELSE
[1210]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
[1216]597                   IF ( PRESENT( ar_2d ) )  THEN
[1210]598
[1392]599                      x_out(0) = CMPLX( ar_2d(0,j), 0.0_wp, KIND=wp )
[1216]600                      DO  i = 1, (nx+1)/2 - 1
[1392]601                         x_out(i) = CMPLX( ar_2d(i,j), ar_2d(nx+1-i,j),        &
602                                           KIND=wp )
[1216]603                      ENDDO
[1392]604                      x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0_wp,      &
605                                               KIND=wp )
[1216]606
607                   ELSE
608
[1392]609                      x_out(0) = CMPLX( ar(0,j,k), 0.0_wp, KIND=wp )
[1216]610                      DO  i = 1, (nx+1)/2 - 1
[1392]611                         x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k), KIND=wp )
[1216]612                      ENDDO
[1392]613                      x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp,       &
614                                               KIND=wp )
[1216]615
616                   ENDIF
617
[1210]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
[1216]625          ENDIF
[1210]626#endif
627
[1106]628       ELSEIF ( fft_method == 'system-specific' )  THEN
629
[1815]630#if defined( __ibm )
[1106]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
[1320]638                   CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1,   &
639                               nau1, aux2, nau2 )
[1106]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
[1342]665                   work(1) = 0.0_wp
666                   work(nx+2) = 0.0_wp
[1106]667
[1320]668                   CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx,      & 
669                               aux3, nau1, aux4, nau2 )
[1106]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
[1342]718                   work(1) = 0.0_wp
719                   work(nx+2) = 0.0_wp
[1106]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
[1682]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 
[1106]746    SUBROUTINE fft_x_1d( ar, direction )
747
748
749       IMPLICIT NONE
750
[1682]751       CHARACTER (LEN=*) ::  direction  !<
[1320]752       
[1682]753       INTEGER(iwp) ::  i               !<
754       INTEGER(iwp) ::  ishape(1)       !<
[1]755
[1682]756       LOGICAL ::  forward_fft          !<
[1106]757
[1682]758       REAL(wp), DIMENSION(0:nx)   ::  ar     !<
759       REAL(wp), DIMENSION(0:nx+2) ::  work   !<
760       REAL(wp), DIMENSION(nx+2)   ::  work1  !<
[1320]761       
[1682]762       COMPLEX(wp), DIMENSION(:), ALLOCATABLE ::  cwork  !<
[1320]763       
[1]764#if defined( __ibm )
[1682]765       REAL(wp), DIMENSION(nau2) ::  aux2       !<
766       REAL(wp), DIMENSION(nau2) ::  aux4       !<
[1]767#elif defined( __nec )
[1682]768       REAL(wp), DIMENSION(6*(nx+1)) ::  work2  !<
[1]769#endif
770
[1106]771       IF ( direction == 'forward' )  THEN
772          forward_fft = .TRUE.
773       ELSE
774          forward_fft = .FALSE.
775       ENDIF
776
[1]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     
[1106]784          IF ( forward_fft )   then
[1]785
786             DO  i = 0, nx
[1392]787                cwork(i) = CMPLX( ar(i), KIND=wp )
[1]788             ENDDO
789             ishape = SHAPE( cwork )
790             CALL FFTN( cwork, ishape )
791             DO  i = 0, (nx+1)/2
[1322]792                ar(i) = REAL( cwork(i), KIND=wp )
[1]793             ENDDO
794             DO  i = 1, (nx+1)/2 - 1
795                ar(nx+1-i) = -AIMAG( cwork(i) )
796             ENDDO
797
798          ELSE
799
[1392]800             cwork(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
[1]801             DO  i = 1, (nx+1)/2 - 1
[1392]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 )
[1]804             ENDDO
[1392]805             cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp, KIND=wp )
[1]806
807             ishape = SHAPE( cwork )
808             CALL FFTN( cwork, ishape, inv = .TRUE. )
809
810             DO  i = 0, nx
[1322]811                ar(i) = REAL( cwork(i), KIND=wp )
[1]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
[1106]823          IF ( forward_fft )  THEN
[1]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
[1342]843             work(1)    = 0.0_wp
844             work(nx+2) = 0.0_wp
[1]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
[1216]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
[1322]860                ar(i) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
[1216]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
[1392]868             x_out(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
[1216]869             DO  i = 1, (nx+1)/2 - 1
[1392]870                x_out(i) = CMPLX( ar(i), ar(nx+1-i), KIND=wp )
[1216]871             ENDDO
[1392]872             x_out((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp, KIND=wp )
[1216]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
[1]880       ELSEIF ( fft_method == 'system-specific' )  THEN
881
[1815]882#if defined( __ibm )
[1106]883          IF ( forward_fft )  THEN
[1]884
[1320]885             CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1,   &
[1]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
[1342]903             work(1) = 0.0_wp
904             work(nx+2) = 0.0_wp
[1]905
[1106]906             CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
[1]907                         aux4, nau2 )
908
909             DO  i = 0, nx
910                ar(i) = work(i)
911             ENDDO
912
913          ENDIF
914#elif defined( __nec )
[1106]915          IF ( forward_fft )  THEN
[1]916
917             work(0:nx) = ar(0:nx)
918
[1106]919             CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
920     
[1]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
[1342]936             work(1) = 0.0_wp
937             work(nx+2) = 0.0_wp
[1]938
[1106]939             CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
[1]940
941             ar(0:nx) = work(0:nx)
942
943          ENDIF
944#endif
945
946       ENDIF
947
[1106]948    END SUBROUTINE fft_x_1d
[1]949
[1682]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 
[1216]974    SUBROUTINE fft_y( ar, direction, ar_tr, nxl_y_bound, nxr_y_bound, nxl_y_l, &
975                      nxr_y_l )
[1]976
977
978       IMPLICIT NONE
979
[1682]980       CHARACTER (LEN=*) ::  direction  !<
[1320]981       
[1682]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      !<
[1106]990
[1682]991       LOGICAL ::  forward_fft  !<
[1106]992
[1682]993       REAL(wp), DIMENSION(0:ny+2) ::  work   !<
994       REAL(wp), DIMENSION(ny+2)   ::  work1  !<
[1320]995       
[1682]996       COMPLEX(wp), DIMENSION(:), ALLOCATABLE ::  cwork  !<
[1320]997       
[1106]998#if defined( __ibm )
[1682]999       REAL(wp), DIMENSION(nau2) ::  auy2  !<
1000       REAL(wp), DIMENSION(nau2) ::  auy4  !<
[1106]1001#elif defined( __nec )
[1682]1002       REAL(wp), DIMENSION(6*(ny+1)) ::  work2  !<
[1106]1003#endif
1004
[1320]1005       REAL(wp), DIMENSION(0:ny,nxl_y_l:nxr_y_l,nzb_y:nzt_y)         ::        &
[1682]1006          ar     !<
[1320]1007       REAL(wp), DIMENSION(0:ny,nxl_y_bound:nxr_y_bound,nzb_y:nzt_y) ::        &
[1682]1008          ar_tr  !<
[1320]1009
[1106]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
[1216]1028                DO  i = nxl_y_l, nxr_y_l
[1106]1029
1030                   DO  j = 0, ny
[1392]1031                      cwork(j) = CMPLX( ar(j,i,k), KIND=wp )
[1106]1032                   ENDDO
1033
1034                   jshape = SHAPE( cwork )
1035                   CALL FFTN( cwork, jshape )
1036
1037                   DO  j = 0, (ny+1)/2
[1322]1038                      ar_tr(j,i,k) = REAL( cwork(j), KIND=wp )
[1106]1039                   ENDDO
1040                   DO  j = 1, (ny+1)/2 - 1
[1216]1041                      ar_tr(ny+1-j,i,k) = -AIMAG( cwork(j) )
[1106]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
[1216]1053                DO  i = nxl_y_l, nxr_y_l
[1106]1054
[1392]1055                   cwork(0) = CMPLX( ar_tr(0,i,k), 0.0_wp, KIND=wp )
[1106]1056                   DO  j = 1, (ny+1)/2 - 1
[1392]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 )
[1106]1061                   ENDDO
[1392]1062                   cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp,       &
1063                                            KIND=wp )
[1106]1064
1065                   jshape = SHAPE( cwork )
1066                   CALL FFTN( cwork, jshape, inv = .TRUE. )
1067
1068                   DO  j = 0, ny
[1322]1069                      ar(j,i,k) = REAL( cwork(j), KIND=wp )
[1106]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
[1304]1087             !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
[1106]1088             !$OMP DO
1089             DO  k = nzb_y, nzt_y
[1216]1090                DO  i = nxl_y_l, nxr_y_l
[1106]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
[1216]1096                      ar_tr(j,i,k) = work(2*j)
[1106]1097                   ENDDO
1098                   DO  j = 1, (ny+1)/2 - 1
[1216]1099                      ar_tr(ny+1-j,i,k) = work(2*j+1)
[1106]1100                   ENDDO
1101
1102                ENDDO
1103             ENDDO
1104             !$OMP END PARALLEL
1105
1106          ELSE
1107
[1304]1108             !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
[1106]1109             !$OMP DO
1110             DO  k = nzb_y, nzt_y
[1216]1111                DO  i = nxl_y_l, nxr_y_l
[1106]1112
1113                   DO  j = 0, (ny+1)/2
[1216]1114                      work(2*j) = ar_tr(j,i,k)
[1106]1115                   ENDDO
1116                   DO  j = 1, (ny+1)/2 - 1
[1216]1117                      work(2*j+1) = ar_tr(ny+1-j,i,k)
[1106]1118                   ENDDO
[1342]1119                   work(1)    = 0.0_wp
1120                   work(ny+2) = 0.0_wp
[1106]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
[1210]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
[1216]1139                DO  i = nxl_y_l, nxr_y_l
[1210]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
[1322]1145                      ar_tr(j,i,k) = REAL( y_out(j), KIND=wp ) / (ny+1)
[1210]1146                   ENDDO
1147                   DO  j = 1, (ny+1)/2 - 1
[1216]1148                      ar_tr(ny+1-j,i,k) = AIMAG( y_out(j) ) / (ny+1)
[1210]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
[1216]1160                DO  i = nxl_y_l, nxr_y_l
[1210]1161
[1392]1162                   y_out(0) = CMPLX( ar_tr(0,i,k), 0.0_wp, KIND=wp )
[1210]1163                   DO  j = 1, (ny+1)/2 - 1
[1398]1164                      y_out(j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k),       &
1165                                        KIND=wp )
[1210]1166                   ENDDO
[1392]1167                   y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp,       &
1168                                            KIND=wp )
[1210]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
[1106]1180       ELSEIF ( fft_method == 'system-specific' )  THEN
1181
[1815]1182#if defined( __ibm )
[1106]1183          IF ( forward_fft)  THEN
1184
1185             !$OMP PARALLEL PRIVATE ( work, i, j, k )
1186             !$OMP DO
1187             DO  k = nzb_y, nzt_y
[1216]1188                DO  i = nxl_y_l, nxr_y_l
[1106]1189
[1320]1190                   CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1,   & 
1191                               nau1, auy2, nau2 )
[1106]1192
1193                   DO  j = 0, (ny+1)/2
[1216]1194                      ar_tr(j,i,k) = work(2*j)
[1106]1195                   ENDDO
1196                   DO  j = 1, (ny+1)/2 - 1
[1216]1197                      ar_tr(ny+1-j,i,k) = work(2*j+1)
[1106]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
[1216]1209                DO  i = nxl_y_l, nxr_y_l
[1106]1210
1211                   DO  j = 0, (ny+1)/2
[1216]1212                      work(2*j) = ar_tr(j,i,k)
[1106]1213                   ENDDO
1214                   DO  j = 1, (ny+1)/2 - 1
[1216]1215                      work(2*j+1) = ar_tr(ny+1-j,i,k)
[1106]1216                   ENDDO
[1342]1217                   work(1)    = 0.0_wp
1218                   work(ny+2) = 0.0_wp
[1106]1219
[1320]1220                   CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny,      &
1221                               auy3, nau1, auy4, nau2 )
[1106]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
[1216]1238                DO  i = nxl_y_l, nxr_y_l
[1106]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
[1216]1245                      ar_tr(j,i,k) = work(2*j)
[1106]1246                   ENDDO
1247                   DO  j = 1, (ny+1)/2 - 1
[1216]1248                      ar_tr(ny+1-j,i,k) = work(2*j+1)
[1106]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
[1216]1260                DO  i = nxl_y_l, nxr_y_l
[1106]1261
1262                   DO  j = 0, (ny+1)/2
[1216]1263                      work(2*j) = ar_tr(j,i,k)
[1106]1264                   ENDDO
1265                   DO  j = 1, (ny+1)/2 - 1
[1216]1266                      work(2*j+1) = ar_tr(ny+1-j,i,k)
[1106]1267                   ENDDO
[1342]1268                   work(1) = 0.0_wp
1269                   work(ny+2) = 0.0_wp
[1106]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
[1682]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 
[1106]1295    SUBROUTINE fft_y_1d( ar, direction )
1296
1297
1298       IMPLICIT NONE
1299
1300       CHARACTER (LEN=*) ::  direction
[1320]1301       
[1682]1302       INTEGER(iwp) ::  j          !<
1303       INTEGER(iwp) ::  jshape(1)  !<
[1]1304
[1682]1305       LOGICAL ::  forward_fft  !<
[1106]1306
[1682]1307       REAL(wp), DIMENSION(0:ny)    ::  ar     !<
1308       REAL(wp), DIMENSION(0:ny+2)  ::  work   !<
1309       REAL(wp), DIMENSION(ny+2)    ::  work1  !<
[1320]1310       
[1682]1311       COMPLEX(wp), DIMENSION(:), ALLOCATABLE ::  cwork  !<
[1320]1312       
[1]1313#if defined( __ibm )
[1682]1314       REAL(wp), DIMENSION(nau2) ::  auy2  !<
1315       REAL(wp), DIMENSION(nau2) ::  auy4  !<
[1]1316#elif defined( __nec )
[1682]1317       REAL(wp), DIMENSION(6*(ny+1)) ::  work2  !<
[1]1318#endif
1319
[1106]1320       IF ( direction == 'forward' )  THEN
1321          forward_fft = .TRUE.
1322       ELSE
1323          forward_fft = .FALSE.
1324       ENDIF
1325
[1]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
[1106]1333          IF ( forward_fft )  THEN
[1]1334
1335             DO  j = 0, ny
[1392]1336                cwork(j) = CMPLX( ar(j), KIND=wp )
[1]1337             ENDDO
1338
1339             jshape = SHAPE( cwork )
1340             CALL FFTN( cwork, jshape )
1341
1342             DO  j = 0, (ny+1)/2
[1322]1343                ar(j) = REAL( cwork(j), KIND=wp )
[1]1344             ENDDO
1345             DO  j = 1, (ny+1)/2 - 1
1346                ar(ny+1-j) = -AIMAG( cwork(j) )
1347             ENDDO
1348
1349          ELSE
1350
[1392]1351             cwork(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
[1]1352             DO  j = 1, (ny+1)/2 - 1
[1392]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 )
[1]1355             ENDDO
[1392]1356             cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp, KIND=wp )
[1]1357
1358             jshape = SHAPE( cwork )
1359             CALL FFTN( cwork, jshape, inv = .TRUE. )
1360
1361             DO  j = 0, ny
[1322]1362                ar(j) = REAL( cwork(j), KIND=wp )
[1]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
[1106]1374          IF ( forward_fft )  THEN
[1]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
[1342]1394             work(1)    = 0.0_wp
1395             work(ny+2) = 0.0_wp
[1]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
[1216]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
[1322]1411                ar(j) = REAL( y_out(j), KIND=wp ) / (ny+1)
[1216]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
[1392]1419             y_out(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
[1216]1420             DO  j = 1, (ny+1)/2 - 1
[1392]1421                y_out(j) = CMPLX( ar(j), ar(ny+1-j), KIND=wp )
[1216]1422             ENDDO
[1392]1423             y_out((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp, KIND=wp )
[1216]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
[1]1431       ELSEIF ( fft_method == 'system-specific' )  THEN
1432
[1815]1433#if defined( __ibm )
[1106]1434          IF ( forward_fft )  THEN
[1]1435
[1320]1436             CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1,   &
[1]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
[1342]1454             work(1)    = 0.0_wp
1455             work(ny+2) = 0.0_wp
[1]1456
[1320]1457             CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3,      &
1458                         nau1, auy4, nau2 )
[1]1459
1460             DO  j = 0, ny
1461                ar(j) = work(j)
1462             ENDDO
1463
1464          ENDIF
1465#elif defined( __nec )
[1106]1466          IF ( forward_fft )  THEN
[1]1467
1468             work(0:ny) = ar(0:ny)
1469
[1106]1470             CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
[1]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
[1342]1487             work(1) = 0.0_wp
1488             work(ny+2) = 0.0_wp
[1]1489
[1106]1490             CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
[1]1491
1492             ar(0:ny) = work(0:ny)
1493
1494          ENDIF
1495#endif
1496
1497       ENDIF
1498
[1106]1499    END SUBROUTINE fft_y_1d
[1]1500
[1682]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 
[1]1510    SUBROUTINE fft_x_m( ar, direction )
1511
1512
1513       IMPLICIT NONE
1514
[1682]1515       CHARACTER (LEN=*) ::  direction  !<
[1320]1516       
[1682]1517       INTEGER(iwp) ::  i     !<
1518       INTEGER(iwp) ::  k     !<
1519       INTEGER(iwp) ::  siza  !<
[3241]1520#if defined( __nec )
1521       INTEGER(iwp) ::  sizw
1522#endif
[1]1523
[1682]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  !<
[1320]1527       
[3241]1528#if defined( __nec )
1529       COMPLEX(wp), DIMENSION(:,:), ALLOCATABLE ::  work
1530#endif
[1]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)
[1342]1539             ai(nx+1:,:)   = 0.0_wp
[1]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
[1342]1561                ai(1,k) = 0.0_wp
1562                ai(nx+2,k) = 0.0_wp
[1]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 )
[2300]1574          ALLOCATE( work((nx+4)/2+1,nz+1) )
[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
[1320]1583             CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4,       &
[1]1584                          trig_xf, work1, 0 )
1585
1586             ai(0:nx,1:nz) = ar(0:nx,1:nz)
1587             IF ( nz1 > nz )  THEN
[1342]1588                ai(:,nz1) = 0.0_wp
[1]1589             ENDIF
1590
[1320]1591             CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw,         &
[1]1592                          trig_xf, work1, 0 )
1593
1594             DO  k = 1, nz
1595                DO  i = 0, (nx+1)/2
[1322]1596                   ar(i,k) = REAL( work(i+1,k), KIND=wp )
[1]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
[1320]1608             CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4,       &
[1]1609                          trig_xb, work1, 0 )
1610
1611             IF ( nz1 > nz )  THEN
[1342]1612                work(:,nz1) = 0.0_wp
[1]1613             ENDIF
1614             DO  k = 1, nz
[1392]1615                work(1,k) = CMPLX( ar(0,k), 0.0_wp, KIND=wp )
[1]1616                DO  i = 1, (nx+1)/2 - 1
[1392]1617                   work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k), KIND=wp )
[1]1618                ENDDO
[1392]1619                work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0_wp, KIND=wp )
[1]1620             ENDDO
1621
[1106]1622             CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, &
[1]1623                          trig_xb, work1, 0 )
1624
1625             ar(0:nx,1:nz) = ai(0:nx,1:nz)
1626
1627          ENDIF
1628
[2300]1629          DEALLOCATE( work )
[1]1630#endif
1631
1632       ENDIF
1633
1634    END SUBROUTINE fft_x_m
1635
[1682]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 
[1]1645    SUBROUTINE fft_y_m( ar, ny1, direction )
1646
1647
1648       IMPLICIT NONE
1649
[1682]1650       CHARACTER (LEN=*) ::  direction  !<
[1320]1651       
[1682]1652       INTEGER(iwp) ::  j     !<
1653       INTEGER(iwp) ::  k     !<
1654       INTEGER(iwp) ::  ny1   !<
1655       INTEGER(iwp) ::  siza  !<
[3241]1656#if defined( __nec )
1657       INTEGER(iwp) ::  sizw
1658#endif
[1]1659
[1682]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  !<
[1]1663
[3241]1664#if defined( __nec )
1665       COMPLEX(wp), DIMENSION(:,:), ALLOCATABLE ::  work
1666#endif
[2300]1667
[3241]1668
[1]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)
[1342]1676             ai(ny+1:,:)   = 0.0_wp
[1]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
[1342]1698                ai(1,k) = 0.0_wp
1699                ai(ny+2,k) = 0.0_wp
[1]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 )
[2300]1711          ALLOCATE( work((ny+4)/2+1,nz+1) )
[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
[1106]1720             CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
[1]1721                          trig_yf, work1, 0 )
1722
1723             ai(0:ny,1:nz) = ar(0:ny,1:nz)
1724             IF ( nz1 > nz )  THEN
[1342]1725                ai(:,nz1) = 0.0_wp
[1]1726             ENDIF
1727
[1106]1728             CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, &
[1]1729                          trig_yf, work1, 0 )
1730
1731             DO  k = 1, nz
1732                DO  j = 0, (ny+1)/2
[1322]1733                   ar(j,k) = REAL( work(j+1,k), KIND=wp )
[1]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
[1106]1745             CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
[1]1746                          trig_yb, work1, 0 )
1747
1748             IF ( nz1 > nz )  THEN
[1342]1749                work(:,nz1) = 0.0_wp
[1]1750             ENDIF
1751             DO  k = 1, nz
[1392]1752                work(1,k) = CMPLX( ar(0,k), 0.0_wp, KIND=wp )
[1]1753                DO  j = 1, (ny+1)/2 - 1
[1392]1754                   work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k), KIND=wp )
[1]1755                ENDDO
[1392]1756                work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0_wp, KIND=wp )
[1]1757             ENDDO
1758
[1106]1759             CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, &
[1]1760                          trig_yb, work1, 0 )
1761
1762             ar(0:ny,1:nz) = ai(0:ny,1:nz)
1763
1764          ENDIF
1765
[2300]1766          DEALLOCATE( work )
[1]1767#endif
1768
1769       ENDIF
1770
1771    END SUBROUTINE fft_y_m
1772
[1106]1773
[1]1774 END MODULE fft_xy
Note: See TracBrowser for help on using the repository browser.