Home

Context Navigation

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

Last change on this file since 1850 was 1850, checked in by maronga, 8 years ago
added _mod string to several filenames to meet the naming convection for modules
Property svn:keywords set to `Id`
File size: 57.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |