Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |