Home

Context Navigation

fft_xy.f90 @ 1402

Last change on this file since 1402 was 1402, checked in by raasch, 11 years ago

output of location messages complemented, output of location bar added
(Makefile, check_parameters, cpulog, init_pegrid, init_3d_model, message, palm, parin, time_integration, new: progress_bar)
preprocessor switch intel_compiler added, -r8 compiler options removed
(.mrun.config.default, .mrun.config.imuk, .mrun.config.kiaps)

batch_job added to envpar-NAMELIST
(mrun, parin)

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: