Home

Context Navigation

source: palm/trunk/SOURCE/fft_xy.f90 @ 1

Last change on this file since 1 was 1, checked in by raasch, 17 years ago
Initial repository layout and content
File size: 24.8 KB

Rev	Line
[1]	1	MODULE fft_xy
	2
	3	!------------------------------------------------------------------------------!
	4	! Actual revisions:
	5	! -----------------
	6	!
	7	!
	8	! Former revisions:
	9	! -----------------
	10	! $Log: fft_xy.f90,v $
	11	! Revision 1.4 2006/03/28 12:27:09 raasch
	12	! Stop when system-specific fft is selected on NEC. For unknown reasons this
	13	! causes a program abort during first allocation in init_grid.
	14	!
	15	! Revision 1.3 2006/02/23 12:22:41 raasch
	16	! System specfic routines are not allowed to be used on ibmy
	17	!
	18	! Revision 1.2 2004/04/30 11:44:27 raasch
	19	! Module renamed from fft_for_1d_decomp to fft_xy, 1d-routines renamed to
	20	! fft_x and fft_y,
	21	! function FFT replaced by subroutine FFTN due to problems with 64-bit
	22	! mode on ibm,
	23	! shape of array cwork is explicitly stored in ishape/jshape and handled
	24	! to routine FFTN instead of shape-function (due to compiler error on
	25	! decalpha),
	26	! non vectorized FFT for nec included
	27	!
	28	! Revision 1.1 2004/04/30 11:37:14 raasch
	29	! Initial revision
	30	!
	31	! Revision 1.6 2003/08/01 08:15:09 raasch
	32	! No abort in fft_1dd_init on t3e-systems if system specific routines are used
	33	!
	34	! Revision 1.5 2003/04/16 12:51:47 raasch
	35	! Temperton-fft implemented in routines fft_x_1dd and fft_y_1dd
	36	!
	37	! Revision 1.4 2003/03/16 09:37:26 raasch
	38	! Two underscores (_) are placed in front of all define-strings
	39	!
	40	! Revision 1.3 2003/03/12 16:29:39 raasch
	41	! Routines fft_x_1dd_m and fft_y_1dd_m added (suitable for multiple 1d-fft on
	42	! vector processors)
	43	!
	44	! Revision 1.2 2002/12/19 14:49:18 raasch
	45	! STOP statement replaced by call of subroutine local_stop
	46	!
	47	! Revision 1.1 2002/06/11 13:00:49 raasch
	48	! Initial revision
	49	!
	50	!
	51	! Description:
	52	! ------------
	53	! Fast Fourier transformation along x and y for 1d domain decomposition along x.
	54	! Original version: Klaus Ketelsen (May 2002)
	55	!------------------------------------------------------------------------------!
	56
	57	USE array_kind
	58	USE control_parameters
	59	USE indices
	60	USE singleton
	61	USE temperton_fft
	62
	63	IMPLICIT NONE
	64
	65	PRIVATE
	66	PUBLIC fft_x, fft_y, fft_init, fft_x_m, fft_y_m
	67
	68	INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: ifax_x, ifax_y
	69
	70	LOGICAL, SAVE :: init_fft = .FALSE.
	71
	72	REAL, SAVE :: sqr_nx, sqr_ny
	73	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trigs_x, trigs_y
	74
	75	#if defined( __ibm )
	76	INTEGER, PARAMETER :: nau1 = 20000, nau2 = 22000
	77	!
	78	!-- The following working arrays contain tables and have to be "save" and
	79	!-- shared in OpenMP sense
	80	REAL, DIMENSION(nau1), SAVE :: aux1, auy1, aux3, auy3
	81	#elif defined( __nec )
	82	INTEGER, SAVE :: nz1
	83	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trig_xb, trig_xf, trig_yb, &
	84	trig_yf
	85	#endif
	86
	87	!
	88	!-- Public interfaces
	89	INTERFACE fft_init
	90	MODULE PROCEDURE fft_init
	91	END INTERFACE fft_init
	92
	93	INTERFACE fft_x
	94	MODULE PROCEDURE fft_x
	95	END INTERFACE fft_x
	96
	97	INTERFACE fft_y
	98	MODULE PROCEDURE fft_y
	99	END INTERFACE fft_y
	100
	101	INTERFACE fft_x_m
	102	MODULE PROCEDURE fft_x_m
	103	END INTERFACE fft_x_m
	104
	105	INTERFACE fft_y_m
	106	MODULE PROCEDURE fft_y_m
	107	END INTERFACE fft_y_m
	108
	109	CONTAINS
	110
	111
	112	SUBROUTINE fft_init
	113
	114	IMPLICIT NONE
	115
	116	!
	117	!-- The following temporary working arrays have to be on stack or private
	118	!-- in OpenMP sense
	119	#if defined( __ibm )
	120	REAL, DIMENSION(0:nx+2) :: workx
	121	REAL, DIMENSION(0:ny+2) :: worky
	122	REAL, DIMENSION(nau2) :: aux2, auy2, aux4, auy4
	123	#elif defined( __nec )
	124	REAL, DIMENSION(0:nx+3,nz+1) :: work_x
	125	REAL, DIMENSION(0:ny+3,nz+1) :: work_y
	126	REAL, DIMENSION(6*(nx+3),nz+1) :: workx
	127	REAL, DIMENSION(6*(ny+3),nz+1) :: worky
	128	#endif
	129
	130	!
	131	!-- Return, if already called
	132	IF ( init_fft ) THEN
	133	RETURN
	134	ELSE
	135	init_fft = .TRUE.
	136	ENDIF
	137
	138	IF ( fft_method == 'system-specific' ) THEN
	139
	140	sqr_nx = SQRT( 1.0 / ( nx + 1.0 ) )
	141	sqr_ny = SQRT( 1.0 / ( ny + 1.0 ) )
	142	#if defined( __ibm ) && ! defined( __ibmy_special )
	143	!
	144	!-- Initialize tables for fft along x
	145	CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1, 1, sqr_nx, aux1, nau1, &
	146	aux2, nau2 )
	147	CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_nx, aux3, nau1, &
	148	aux4, nau2 )
	149	!
	150	!-- Initialize tables for fft along y
	151	CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1, 1, sqr_ny, auy1, nau1, &
	152	auy2, nau2 )
	153	CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_ny, auy3, nau1, &
	154	auy4, nau2 )
	155	#elif defined( __nec )
	156	PRINT*, '+++ fft_init: fft method "', fft_method, &
	157	'" currently does not work on NEC'
	158	CALL local_stop
	159
	160	ALLOCATE( trig_xb(2(nx+1)), trig_xf(2(nx+1)), &
	161	trig_yb(2(ny+1)), trig_yf(2(ny+1)) )
	162
	163	work_x = 0.0
	164	work_y = 0.0
	165	nz1 = nz + MOD( nz+1, 2 ) ! odd nz slows down fft significantly
	166	! when using the NEC ffts
	167
	168	!
	169	!-- Initialize tables for fft along x (non-vector and vector case (M))
	170	CALL DZFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xf, workx, 0 )
	171	CALL ZDFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xb, workx, 0 )
	172	CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
	173	trig_xf, workx, 0 )
	174	CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
	175	trig_xb, workx, 0 )
	176	!
	177	!-- Initialize tables for fft along y (non-vector and vector case (M))
	178	CALL DZFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yf, worky, 0 )
	179	CALL ZDFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yb, worky, 0 )
	180	CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
	181	trig_yf, worky, 0 )
	182	CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
	183	trig_yb, worky, 0 )
	184	#else
	185	PRINT*, '+++ fft_init: no system-specific fft-call available'
	186	CALL local_stop
	187	#endif
	188	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	189	!
	190	!-- Temperton-algorithm
	191	!-- Initialize tables for fft along x and y
	192	ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
	193
	194	CALL set99( trigs_x, ifax_x, nx+1 )
	195	CALL set99( trigs_y, ifax_y, ny+1 )
	196
	197	ELSEIF ( fft_method == 'singleton-algorithm' ) THEN
	198
	199	CONTINUE
	200
	201	ELSE
	202
	203	PRINT*, '+++ fft_init: fft method "', fft_method, &
	204	'" not available'
	205	CALL local_stop
	206
	207	ENDIF
	208
	209	END SUBROUTINE fft_init
	210
	211
	212	SUBROUTINE fft_x( ar, direction )
	213
	214	!----------------------------------------------------------------------!
	215	! fft_x !
	216	! !
	217	! Fourier-transformation along x-direction !
	218	! !
	219	! fft_x uses internal algorithms (Singleton or Temperton) or !
	220	! system-specific routines, if they are available !
	221	!----------------------------------------------------------------------!
	222
	223	IMPLICIT NONE
	224
	225	CHARACTER (LEN=*) :: direction
	226	INTEGER :: i, ishape(1)
	227
	228	!kk REAL, DIMENSION(:) :: ar !kk Does NOT work (Bug??)
	229	REAL, DIMENSION(0:nx) :: ar
	230	REAL, DIMENSION(0:nx+2) :: work
	231	REAL, DIMENSION(nx+2) :: work1
	232	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
	233	#if defined( __ibm )
	234	REAL, DIMENSION(nau2) :: aux2, aux4
	235	#elif defined( __nec )
	236	REAL, DIMENSION(6*(nx+1)) :: work2
	237	#endif
	238
	239	IF ( fft_method == 'singleton-algorithm' ) THEN
	240
	241	!
	242	!-- Performing the fft with singleton's software works on every system,
	243	!-- since it is part of the model
	244	ALLOCATE( cwork(0:nx) )
	245
	246	IF ( direction == 'forward') then
	247
	248	DO i = 0, nx
	249	cwork(i) = CMPLX( ar(i) )
	250	ENDDO
	251	ishape = SHAPE( cwork )
	252	CALL FFTN( cwork, ishape )
	253
	254	DO i = 0, (nx+1)/2
	255	ar(i) = REAL( cwork(i) )
	256	ENDDO
	257	DO i = 1, (nx+1)/2 - 1
	258	ar(nx+1-i) = -AIMAG( cwork(i) )
	259	ENDDO
	260
	261	ELSE
	262
	263	cwork(0) = CMPLX( ar(0), 0.0 )
	264	DO i = 1, (nx+1)/2 - 1
	265	cwork(i) = CMPLX( ar(i), -ar(nx+1-i) )
	266	cwork(nx+1-i) = CMPLX( ar(i), ar(nx+1-i) )
	267	ENDDO
	268	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
	269
	270	ishape = SHAPE( cwork )
	271	CALL FFTN( cwork, ishape, inv = .TRUE. )
	272
	273	DO i = 0, nx
	274	ar(i) = REAL( cwork(i) )
	275	ENDDO
	276
	277	ENDIF
	278
	279	DEALLOCATE( cwork )
	280
	281	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	282
	283	!
	284	!-- Performing the fft with Temperton's software works on every system,
	285	!-- since it is part of the model
	286	IF ( direction == 'forward' ) THEN
	287
	288	work(0:nx) = ar
	289	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
	290
	291	DO i = 0, (nx+1)/2
	292	ar(i) = work(2*i)
	293	ENDDO
	294	DO i = 1, (nx+1)/2 - 1
	295	ar(nx+1-i) = work(2*i+1)
	296	ENDDO
	297
	298	ELSE
	299
	300	DO i = 0, (nx+1)/2
	301	work(2*i) = ar(i)
	302	ENDDO
	303	DO i = 1, (nx+1)/2 - 1
	304	work(2*i+1) = ar(nx+1-i)
	305	ENDDO
	306	work(1) = 0.0
	307	work(nx+2) = 0.0
	308
	309	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
	310	ar = work(0:nx)
	311
	312	ENDIF
	313
	314	ELSEIF ( fft_method == 'system-specific' ) THEN
	315
	316	#if defined( __ibm ) && ! defined( __ibmy_special )
	317	IF ( direction == 'forward' ) THEN
	318
	319	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_nx, aux1, nau1, &
	320	aux2, nau2 )
	321
	322	DO i = 0, (nx+1)/2
	323	ar(i) = work(2*i)
	324	ENDDO
	325	DO i = 1, (nx+1)/2 - 1
	326	ar(nx+1-i) = work(2*i+1)
	327	ENDDO
	328
	329	ELSE
	330
	331	DO i = 0, (nx+1)/2
	332	work(2*i) = ar(i)
	333	ENDDO
	334	DO i = 1, (nx+1)/2 - 1
	335	work(2*i+1) = ar(nx+1-i)
	336	ENDDO
	337	work(1) = 0.0
	338	work(nx+2) = 0.0
	339
	340	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_nx, aux3, nau1, &
	341	aux4, nau2 )
	342
	343	DO i = 0, nx
	344	ar(i) = work(i)
	345	ENDDO
	346
	347	ENDIF
	348	#elif defined( __nec )
	349	IF ( direction == 'forward' ) THEN
	350
	351	work(0:nx) = ar(0:nx)
	352
	353	CALL DZFFT( 1, nx+1, sqr_nx, work, work, trig_xf, work2, 0 )
	354
	355	DO i = 0, (nx+1)/2
	356	ar(i) = work(2*i)
	357	ENDDO
	358	DO i = 1, (nx+1)/2 - 1
	359	ar(nx+1-i) = work(2*i+1)
	360	ENDDO
	361
	362	ELSE
	363
	364	DO i = 0, (nx+1)/2
	365	work(2*i) = ar(i)
	366	ENDDO
	367	DO i = 1, (nx+1)/2 - 1
	368	work(2*i+1) = ar(nx+1-i)
	369	ENDDO
	370	work(1) = 0.0
	371	work(nx+2) = 0.0
	372
	373	CALL ZDFFT( -1, nx+1, sqr_nx, work, work, trig_xb, work2, 0 )
	374
	375	ar(0:nx) = work(0:nx)
	376
	377	ENDIF
	378	#else
	379	PRINT*, '+++ fft_x: no system-specific fft-call available'
	380	CALL local_stop
	381	#endif
	382	ELSE
	383
	384	PRINT*, '+++ fft_x: fft method "', fft_method, '" not available'
	385	CALL local_stop
	386
	387	ENDIF
	388
	389	END SUBROUTINE fft_x
	390
	391	SUBROUTINE fft_y( ar, direction )
	392
	393	!----------------------------------------------------------------------!
	394	! fft_y !
	395	! !
	396	! Fourier-transformation along y-direction !
	397	! !
	398	! fft_y uses internal algorithms (Singleton or Temperton) or !
	399	! system-specific routines, if they are available !
	400	!----------------------------------------------------------------------!
	401
	402	IMPLICIT NONE
	403
	404	CHARACTER (LEN=*) :: direction
	405	INTEGER :: j, jshape(1)
	406
	407	!kk REAL, DIMENSION(:) :: ar !kk Does NOT work (Bug??)
	408	REAL, DIMENSION(0:ny) :: ar
	409	REAL, DIMENSION(0:ny+2) :: work
	410	REAL, DIMENSION(ny+2) :: work1
	411	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
	412	#if defined( __ibm )
	413	REAL, DIMENSION(nau2) :: auy2, auy4
	414	#elif defined( __nec )
	415	REAL, DIMENSION(6*(ny+1)) :: work2
	416	#endif
	417
	418	IF ( fft_method == 'singleton-algorithm' ) THEN
	419
	420	!
	421	!-- Performing the fft with singleton's software works on every system,
	422	!-- since it is part of the model
	423	ALLOCATE( cwork(0:ny) )
	424
	425	IF ( direction == 'forward') THEN
	426
	427	DO j = 0, ny
	428	cwork(j) = CMPLX( ar(j) )
	429	ENDDO
	430
	431	jshape = SHAPE( cwork )
	432	CALL FFTN( cwork, jshape )
	433
	434	DO j = 0, (ny+1)/2
	435	ar(j) = REAL( cwork(j) )
	436	ENDDO
	437	DO j = 1, (ny+1)/2 - 1
	438	ar(ny+1-j) = -AIMAG( cwork(j) )
	439	ENDDO
	440
	441	ELSE
	442
	443	cwork(0) = CMPLX( ar(0), 0.0 )
	444	DO j = 1, (ny+1)/2 - 1
	445	cwork(j) = CMPLX( ar(j), -ar(ny+1-j) )
	446	cwork(ny+1-j) = CMPLX( ar(j), ar(ny+1-j) )
	447	ENDDO
	448	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
	449
	450	jshape = SHAPE( cwork )
	451	CALL FFTN( cwork, jshape, inv = .TRUE. )
	452
	453	DO j = 0, ny
	454	ar(j) = REAL( cwork(j) )
	455	ENDDO
	456
	457	ENDIF
	458
	459	DEALLOCATE( cwork )
	460
	461	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	462
	463	!
	464	!-- Performing the fft with Temperton's software works on every system,
	465	!-- since it is part of the model
	466	IF ( direction == 'forward' ) THEN
	467
	468	work(0:ny) = ar
	469	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
	470
	471	DO j = 0, (ny+1)/2
	472	ar(j) = work(2*j)
	473	ENDDO
	474	DO j = 1, (ny+1)/2 - 1
	475	ar(ny+1-j) = work(2*j+1)
	476	ENDDO
	477
	478	ELSE
	479
	480	DO j = 0, (ny+1)/2
	481	work(2*j) = ar(j)
	482	ENDDO
	483	DO j = 1, (ny+1)/2 - 1
	484	work(2*j+1) = ar(ny+1-j)
	485	ENDDO
	486	work(1) = 0.0
	487	work(ny+2) = 0.0
	488
	489	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
	490	ar = work(0:ny)
	491
	492	ENDIF
	493
	494	ELSEIF ( fft_method == 'system-specific' ) THEN
	495
	496	#if defined( __ibm ) && ! defined( __ibmy_special )
	497	IF ( direction == 'forward') THEN
	498
	499	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_ny, auy1, nau1, &
	500	auy2, nau2 )
	501
	502	DO j = 0, (ny+1)/2
	503	ar(j) = work(2*j)
	504	ENDDO
	505	DO j = 1, (ny+1)/2 - 1
	506	ar(ny+1-j) = work(2*j+1)
	507	ENDDO
	508
	509	ELSE
	510
	511	DO j = 0, (ny+1)/2
	512	work(2*j) = ar(j)
	513	ENDDO
	514	DO j = 1, (ny+1)/2 - 1
	515	work(2*j+1) = ar(ny+1-j)
	516	ENDDO
	517	work(1) = 0.0
	518	work(ny+2) = 0.0
	519
	520	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_ny, auy3, nau1, &
	521	auy4, nau2 )
	522
	523	DO j = 0, ny
	524	ar(j) = work(j)
	525	ENDDO
	526
	527	ENDIF
	528	#elif defined( __nec )
	529	IF ( direction == 'forward' ) THEN
	530
	531	work(0:ny) = ar(0:ny)
	532
	533	CALL DZFFT( 1, ny+1, sqr_ny, work, work, trig_yf, work2, 0 )
	534
	535	DO j = 0, (ny+1)/2
	536	ar(j) = work(2*j)
	537	ENDDO
	538	DO j = 1, (ny+1)/2 - 1
	539	ar(ny+1-j) = work(2*j+1)
	540	ENDDO
	541
	542	ELSE
	543
	544	DO j = 0, (ny+1)/2
	545	work(2*j) = ar(j)
	546	ENDDO
	547	DO j = 1, (ny+1)/2 - 1
	548	work(2*j+1) = ar(ny+1-j)
	549	ENDDO
	550	work(1) = 0.0
	551	work(ny+2) = 0.0
	552
	553	CALL ZDFFT( -1, ny+1, sqr_ny, work, work, trig_yb, work2, 0 )
	554
	555	ar(0:ny) = work(0:ny)
	556
	557	ENDIF
	558	#else
	559	PRINT*, '+++ fft_y: no system-specific fft-call available'
	560	CALL local_stop
	561	#endif
	562
	563	ELSE
	564
	565	PRINT*, '+++ fft_y: fft method "', fft_method, '" not available'
	566	CALL local_stop
	567
	568	ENDIF
	569
	570	END SUBROUTINE fft_y
	571
	572	SUBROUTINE fft_x_m( ar, direction )
	573
	574	!----------------------------------------------------------------------!
	575	! fft_x_m !
	576	! !
	577	! Fourier-transformation along x-direction !
	578	! Version for 1d domain decomposition !
	579	! using multiple 1D FFT from Math Keisan on NEC !
	580	! or Temperton-algorithm !
	581	! (no singleton-algorithm on NEC because it does not vectorize) !
	582	! !
	583	!----------------------------------------------------------------------!
	584
	585	IMPLICIT NONE
	586
	587	CHARACTER (LEN=*) :: direction
	588	INTEGER :: i, k, siza, sizw
	589
	590	REAL, DIMENSION(0:nx,nz) :: ar
	591	REAL, DIMENSION(0:nx+3,nz+1) :: ai
	592	REAL, DIMENSION(6*(nx+4),nz+1) :: work1
	593	#if defined( __nec )
	594	COMPLEX, DIMENSION((nx+4)/2+1,nz+1) :: work
	595	#endif
	596
	597	IF ( fft_method == 'temperton-algorithm' ) THEN
	598
	599	siza = SIZE( ai, 1 )
	600
	601	IF ( direction == 'forward') THEN
	602
	603	ai(0:nx,1:nz) = ar(0:nx,1:nz)
	604	ai(nx+1:,:) = 0.0
	605
	606	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
	607
	608	DO k = 1, nz
	609	DO i = 0, (nx+1)/2
	610	ar(i,k) = ai(2*i,k)
	611	ENDDO
	612	DO i = 1, (nx+1)/2 - 1
	613	ar(nx+1-i,k) = ai(2*i+1,k)
	614	ENDDO
	615	ENDDO
	616
	617	ELSE
	618
	619	DO k = 1, nz
	620	DO i = 0, (nx+1)/2
	621	ai(2*i,k) = ar(i,k)
	622	ENDDO
	623	DO i = 1, (nx+1)/2 - 1
	624	ai(2*i+1,k) = ar(nx+1-i,k)
	625	ENDDO
	626	ai(1,k) = 0.0
	627	ai(nx+2,k) = 0.0
	628	ENDDO
	629
	630	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
	631
	632	ar(0:nx,1:nz) = ai(0:nx,1:nz)
	633
	634	ENDIF
	635
	636	ELSEIF ( fft_method == 'system-specific' ) THEN
	637
	638	#if defined( __nec )
	639	siza = SIZE( ai, 1 )
	640	sizw = SIZE( work, 1 )
	641
	642	IF ( direction == 'forward') THEN
	643
	644	!
	645	!-- Tables are initialized once more. This call should not be
	646	!-- necessary, but otherwise program aborts in asymmetric case
	647	CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
	648	trig_xf, work1, 0 )
	649
	650	ai(0:nx,1:nz) = ar(0:nx,1:nz)
	651	IF ( nz1 > nz ) THEN
	652	ai(:,nz1) = 0.0
	653	ENDIF
	654
	655	CALL DZFFTM( 1, nx+1, nz1, sqr_nx, ai, siza, work, sizw, &
	656	trig_xf, work1, 0 )
	657
	658	DO k = 1, nz
	659	DO i = 0, (nx+1)/2
	660	ar(i,k) = REAL( work(i+1,k) )
	661	ENDDO
	662	DO i = 1, (nx+1)/2 - 1
	663	ar(nx+1-i,k) = AIMAG( work(i+1,k) )
	664	ENDDO
	665	ENDDO
	666
	667	ELSE
	668
	669	!
	670	!-- Tables are initialized once more. This call should not be
	671	!-- necessary, but otherwise program aborts in asymmetric case
	672	CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
	673	trig_xb, work1, 0 )
	674
	675	IF ( nz1 > nz ) THEN
	676	work(:,nz1) = 0.0
	677	ENDIF
	678	DO k = 1, nz
	679	work(1,k) = CMPLX( ar(0,k), 0.0 )
	680	DO i = 1, (nx+1)/2 - 1
	681	work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
	682	ENDDO
	683	work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
	684	ENDDO
	685
	686	CALL ZDFFTM( -1, nx+1, nz1, sqr_nx, work, sizw, ai, siza, &
	687	trig_xb, work1, 0 )
	688
	689	ar(0:nx,1:nz) = ai(0:nx,1:nz)
	690
	691	ENDIF
	692
	693	#else
	694	PRINT*, '+++ fft_x_m: no system-specific fft-call available'
	695	STOP
	696	#endif
	697
	698	ELSE
	699
	700	PRINT*, '+++ fft_x_m: fft method "', fft_method, '" not available'
	701	CALL local_stop
	702
	703	ENDIF
	704
	705	END SUBROUTINE fft_x_m
	706
	707	SUBROUTINE fft_y_m( ar, ny1, direction )
	708
	709	!----------------------------------------------------------------------!
	710	! fft_y_m !
	711	! !
	712	! Fourier-transformation along y-direction !
	713	! Version for 1d domain decomposition !
	714	! using multiple 1D FFT from Math Keisan on NEC !
	715	! or Temperton-algorithm !
	716	! (no singleton-algorithm on NEC because it does not vectorize) !
	717	! !
	718	!----------------------------------------------------------------------!
	719
	720	IMPLICIT NONE
	721
	722	CHARACTER (LEN=*) :: direction
	723	INTEGER :: j, k, ny1, siza, sizw
	724
	725	REAL, DIMENSION(0:ny1,nz) :: ar
	726	REAL, DIMENSION(0:ny+3,nz+1) :: ai
	727	REAL, DIMENSION(6*(ny+4),nz+1) :: work1
	728	#if defined( __nec )
	729	COMPLEX, DIMENSION((ny+4)/2+1,nz+1) :: work
	730	#endif
	731
	732	IF ( fft_method == 'temperton-algorithm' ) THEN
	733
	734	siza = SIZE( ai, 1 )
	735
	736	IF ( direction == 'forward') THEN
	737
	738	ai(0:ny,1:nz) = ar(0:ny,1:nz)
	739	ai(ny+1:,:) = 0.0
	740
	741	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
	742
	743	DO k = 1, nz
	744	DO j = 0, (ny+1)/2
	745	ar(j,k) = ai(2*j,k)
	746	ENDDO
	747	DO j = 1, (ny+1)/2 - 1
	748	ar(ny+1-j,k) = ai(2*j+1,k)
	749	ENDDO
	750	ENDDO
	751
	752	ELSE
	753
	754	DO k = 1, nz
	755	DO j = 0, (ny+1)/2
	756	ai(2*j,k) = ar(j,k)
	757	ENDDO
	758	DO j = 1, (ny+1)/2 - 1
	759	ai(2*j+1,k) = ar(ny+1-j,k)
	760	ENDDO
	761	ai(1,k) = 0.0
	762	ai(ny+2,k) = 0.0
	763	ENDDO
	764
	765	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
	766
	767	ar(0:ny,1:nz) = ai(0:ny,1:nz)
	768
	769	ENDIF
	770
	771	ELSEIF ( fft_method == 'system-specific' ) THEN
	772
	773	#if defined( __nec )
	774	siza = SIZE( ai, 1 )
	775	sizw = SIZE( work, 1 )
	776
	777	IF ( direction == 'forward') THEN
	778
	779	!
	780	!-- Tables are initialized once more. This call should not be
	781	!-- necessary, but otherwise program aborts in asymmetric case
	782	CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
	783	trig_yf, work1, 0 )
	784
	785	ai(0:ny,1:nz) = ar(0:ny,1:nz)
	786	IF ( nz1 > nz ) THEN
	787	ai(:,nz1) = 0.0
	788	ENDIF
	789
	790	CALL DZFFTM( 1, ny+1, nz1, sqr_ny, ai, siza, work, sizw, &
	791	trig_yf, work1, 0 )
	792
	793	DO k = 1, nz
	794	DO j = 0, (ny+1)/2
	795	ar(j,k) = REAL( work(j+1,k) )
	796	ENDDO
	797	DO j = 1, (ny+1)/2 - 1
	798	ar(ny+1-j,k) = AIMAG( work(j+1,k) )
	799	ENDDO
	800	ENDDO
	801
	802	ELSE
	803
	804	!
	805	!-- Tables are initialized once more. This call should not be
	806	!-- necessary, but otherwise program aborts in asymmetric case
	807	CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
	808	trig_yb, work1, 0 )
	809
	810	IF ( nz1 > nz ) THEN
	811	work(:,nz1) = 0.0
	812	ENDIF
	813	DO k = 1, nz
	814	work(1,k) = CMPLX( ar(0,k), 0.0 )
	815	DO j = 1, (ny+1)/2 - 1
	816	work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
	817	ENDDO
	818	work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
	819	ENDDO
	820
	821	CALL ZDFFTM( -1, ny+1, nz1, sqr_ny, work, sizw, ai, siza, &
	822	trig_yb, work1, 0 )
	823
	824	ar(0:ny,1:nz) = ai(0:ny,1:nz)
	825
	826	ENDIF
	827
	828	#else
	829	PRINT*, '+++ fft_y_m: no system-specific fft-call available'
	830	STOP
	831	#endif
	832
	833	ELSE
	834
	835	PRINT*, '+++ fft_y_m: fft method "', fft_method, '" not available'
	836	CALL local_stop
	837
	838	ENDIF
	839
	840	END SUBROUTINE fft_y_m
	841
	842	END MODULE fft_xy

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |