Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |