1 | |
---|
2 | ! $Id: lmdz_fft_fftw.F90 5134 2024-07-26 15:56:37Z fairhead $ |
---|
3 | |
---|
4 | MODULE lmdz_fft_fftw |
---|
5 | |
---|
6 | #ifdef FFT_FFTW |
---|
7 | |
---|
8 | REAL, SAVE :: scale_factor |
---|
9 | INTEGER, SAVE :: vsize |
---|
10 | INTEGER, PARAMETER :: inc=1 |
---|
11 | |
---|
12 | INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_forward |
---|
13 | INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_backward |
---|
14 | |
---|
15 | CONTAINS |
---|
16 | |
---|
17 | SUBROUTINE Init_fft(iim,nvectmax) |
---|
18 | IMPLICIT NONE |
---|
19 | #INCLUDE <fftw3.f> |
---|
20 | INTEGER :: iim |
---|
21 | INTEGER :: nvectmax |
---|
22 | |
---|
23 | INTEGER :: itmp |
---|
24 | |
---|
25 | INTEGER :: rank |
---|
26 | INTEGER :: howmany |
---|
27 | INTEGER :: istride, idist |
---|
28 | INTEGER :: ostride, odist |
---|
29 | INTEGER, DIMENSION(1) :: n_array, inembed, onembed |
---|
30 | |
---|
31 | REAL, DIMENSION(iim+1,nvectmax) :: dbidon |
---|
32 | COMPLEX, DIMENSION(iim/2+1,nvectmax) :: cbidon |
---|
33 | |
---|
34 | vsize = iim |
---|
35 | scale_factor = 1./SQRT(1.*vsize) |
---|
36 | |
---|
37 | dbidon = 0 |
---|
38 | cbidon = 0 |
---|
39 | |
---|
40 | ALLOCATE(plan_forward(nvectmax)) |
---|
41 | ALLOCATE(plan_backward(nvectmax)) |
---|
42 | |
---|
43 | WRITE(*,*)"!---------------------!" |
---|
44 | WRITE(*,*)"! !" |
---|
45 | WRITE(*,*)"! INITIALISATION FFTW !" |
---|
46 | WRITE(*,*)"! !" |
---|
47 | WRITE(*,*)"!---------------------!" |
---|
48 | |
---|
49 | ! On initialise tous les plans |
---|
50 | DO itmp = 1, nvectmax |
---|
51 | rank = 1 |
---|
52 | n_array(1) = iim |
---|
53 | howmany = itmp |
---|
54 | inembed(1) = iim + 1 ; onembed(1) = iim/2 + 1 |
---|
55 | istride = 1 ; ostride = 1 |
---|
56 | idist = iim + 1 ; odist = iim/2 + 1 |
---|
57 | |
---|
58 | CALL dfftw_plan_many_dft_r2c(plan_forward(itmp), rank, n_array, howmany, & |
---|
59 | dbidon, inembed, istride, idist, & |
---|
60 | cbidon, onembed, ostride, odist, & |
---|
61 | FFTW_ESTIMATE) |
---|
62 | |
---|
63 | rank = 1 |
---|
64 | n_array(1) = iim |
---|
65 | howmany = itmp |
---|
66 | inembed(1) = iim/2 + 1 ; onembed(1) = iim + 1 |
---|
67 | istride = 1 ; ostride = 1 |
---|
68 | idist = iim/2 + 1 ; odist = iim + 1 |
---|
69 | CALL dfftw_plan_many_dft_c2r(plan_backward(itmp), rank, n_array, howmany, & |
---|
70 | cbidon, inembed, istride, idist, & |
---|
71 | dbidon, onembed, ostride, odist, & |
---|
72 | FFTW_ESTIMATE) |
---|
73 | |
---|
74 | ENDDO |
---|
75 | |
---|
76 | WRITE(*,*)"!-------------------------!" |
---|
77 | WRITE(*,*)"! !" |
---|
78 | WRITE(*,*)"! FIN INITIALISATION FFTW !" |
---|
79 | WRITE(*,*)"! !" |
---|
80 | WRITE(*,*)"!-------------------------!" |
---|
81 | |
---|
82 | END SUBROUTINE Init_fft |
---|
83 | |
---|
84 | |
---|
85 | SUBROUTINE fft_forward(vect,TF_vect,nb_vect) |
---|
86 | IMPLICIT NONE |
---|
87 | #INCLUDE <fftw3.f> |
---|
88 | INTEGER,INTENT(IN) :: nb_vect |
---|
89 | REAL,INTENT(IN) :: vect(vsize+inc,nb_vect) |
---|
90 | COMPLEX,INTENT(OUT) :: TF_vect(vsize/2+1,nb_vect) |
---|
91 | |
---|
92 | CALL dfftw_execute_dft_r2c(plan_forward(nb_vect),vect,TF_vect) |
---|
93 | |
---|
94 | TF_vect = scale_factor * TF_vect |
---|
95 | |
---|
96 | END SUBROUTINE fft_forward |
---|
97 | |
---|
98 | SUBROUTINE fft_backward(TF_vect,vect,nb_vect) |
---|
99 | IMPLICIT NONE |
---|
100 | #INCLUDE <fftw3.f> |
---|
101 | INTEGER,INTENT(IN) :: nb_vect |
---|
102 | REAL,INTENT(OUT) :: vect(vsize+inc,nb_vect) |
---|
103 | COMPLEX,INTENT(IN ) :: TF_vect(vsize/2+1,nb_vect) |
---|
104 | |
---|
105 | CALL dfftw_execute_dft_c2r(plan_backward(nb_vect),TF_vect,vect) |
---|
106 | |
---|
107 | vect = scale_factor * vect |
---|
108 | |
---|
109 | END SUBROUTINE fft_backward |
---|
110 | |
---|
111 | #endif |
---|
112 | |
---|
113 | END MODULE lmdz_fft_fftw |
---|