-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfft.f90
116 lines (66 loc) · 3.71 KB
/
fft.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
MODULE fft
USE parameters
IMPLICIT NONE
integer, parameter :: fftw_forward = -1, fftw_backward = 1
integer, parameter :: fftw_real_to_complex = -1, fftw_complex_to_real = 1
integer, parameter :: fftw_estimate = 0, fftw_measure=1,fftw_use_wisdom = 16
integer, parameter :: fftw_in_place = 8
!FFT plans
integer*8 :: plan_r2c
integer*8 :: plan_c2r
integer, parameter :: norm_rc = n1*n2
!For z-spectra.
integer*8 :: plan_ztokz,plan_kztoz
integer, parameter :: fftw_redft10 = 5, fftw_redft01 = 4
CONTAINS
SUBROUTINE initialize_fftw(array2dr,array2di,fr_even,fk_even)
!2D horizontal transforms
double precision, dimension(n1d,n2d) :: array2dr
double complex, dimension(iktx,ikty) :: array2di
!For z-spectra
double precision, dimension(n3) :: fr_even,fk_even
!2D horizontal transforms
call dfftw_plan_dft_r2c_2d(plan_r2c,n1,n2,array2dr,array2di,fftw_estimate+fftw_in_place)
call dfftw_plan_dft_c2r_2d(plan_c2r,n1,n2,array2di,array2dr,fftw_estimate+fftw_in_place)
!For z-spectra
call dfftw_plan_r2r_1d(plan_ztokz,n3,fr_even,fk_even, FFTW_REDFT10, fftw_estimate+fftw_in_place)
call dfftw_plan_r2r_1d(plan_kztoz,n3,fk_even,fr_even, FFTW_REDFT01, fftw_estimate+fftw_in_place)
END SUBROUTINE initialize_fftw
SUBROUTINE fft_r2c(rvar,kvar,nz)
double precision :: rvar(n1d,n2d,nz)
double complex :: kvar(iktx,ikty,nz)
integer :: nz
do iz=1,nz
call dfftw_execute_dft_r2c(plan_r2c,rvar(1,1,iz),kvar(1,1,iz))
end do
kvar=kvar/norm_rc
END SUBROUTINE fft_r2c
SUBROUTINE fft_c2r(kvar,rvar,nz)
double precision :: rvar(n1d,n2d,nz)
double complex :: kvar(iktx,ikty,nz)
integer :: nz
if(dealiasing==1) call realit(kvar,nz)
do iz=1,nz
call dfftw_execute_dft_c2r(plan_c2r,kvar(1,1,iz),rvar(1,1,iz))
end do
END SUBROUTINE fft_c2r
SUBROUTINE realit(kvar,nz) !This routine adds the right modes (satisfying the reality condition) in (0,ky-) so that the r-field is the right one.
double complex, dimension(iktx,ikty,nz) :: kvar
integer :: nz
integer :: iz !dummy index
integer :: inky !index for negative ky's
!Modes (0,-ky) are just the conjugate of (0,ky+) modes
!Recall kya=[ 0 1 2 ... N/2-1 0 -N/2+1 -N/2+2 ... -2 -1]
do iz=1,nz
do iky=2,ikty/2 !n2 must be even, otherwise it won't work.
inky=ikty-iky+2
kvar(1,inky,iz)=dble(kvar(1,iky,iz)) - i*dimag(kvar(1,iky,iz)) ! u(0,ky-)=u*(0,ky+)
end do
end do
END SUBROUTINE realit
SUBROUTINE kill_fftw
!*** r-space <=> h-space ***!
call dfftw_destroy_plan(plan_r2c)
call dfftw_destroy_plan(plan_c2r)
END SUBROUTINE kill_fftw
END MODULE fft