Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |