Home

Context Navigation

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

Last change on this file since 3013 was 2718, checked in by maronga, 7 years ago
deleting of deprecated files; headers updated where needed
Property svn:keywords set to `Id`
File size: 53.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |