Home

Context Navigation

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

Last change on this file since 4649 was 4646, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 59.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |