Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |