Home

Context Navigation

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

Last change on this file since 1724 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 57.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |