Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |