Home

Context Navigation

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

Last change on this file since 4180 was 4180, checked in by scharf, 5 years ago
removed comments in 'Former revisions' section that are older than 01.01.2019
Property svn:keywords set to `Id`
File size: 52.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |