Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |