Home

Context Navigation

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

Last change on this file since 4174 was 4069, checked in by Giersch, 5 years ago
Bugfix for masked output, compiler warning removed, test case for wind turbine model revised
Property svn:keywords set to `Id`
File size: 55.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |