Home

Context Navigation

source: palm/trunk/SOURCE/fft_xy.f90 @ 1375

Last change on this file since 1375 was 1375, checked in by raasch, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 57.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |