Home

Context Navigation

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

Last change on this file since 2198 was 2119, checked in by raasch, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 52.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |