Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |