Home

Context Navigation

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

Last change on this file since 4366 was 4366, checked in by raasch, 4 years ago
code vectorization for NEC Aurora: vectorized version of Temperton FFT, vectorization of Newtor iteration for calculating the Obukhov length
Property svn:keywords set to `Id`
File size: 60.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |