Home

Context Navigation

fft_xy_mod.f90 @ 4830

Last change on this file since 4830 was 4828, checked in by Giersch, 4 years ago
Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code
Property svn:keywords set to `Id`
File size: 53.2 KB

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

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

Context Navigation

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

Download in other formats: