source: LMDZ6/branches/Amaury_dev/libf/filtrez/lmdz_fft_fftw.F90 @ 5112

Last change on this file since 5112 was 5107, checked in by abarral, 4 months ago

Turn coefils.h into lmdz_coefils.f90
Put filtreg.F90, inifgn.f90, jacobi.F90, eigen_sort.f90, acc.f90 inside lmdz_filtreg.F90
Turn mod_* into lmdz_* in filtrez
Delete obsolete parafilt.h*
(lint) remove spaces between routine name and args
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.0 KB
Line 
1
2! $Id: lmdz_fft_fftw.F90 5107 2024-07-24 08:26:10Z abarral $
3
4MODULE 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 
15CONTAINS
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 
113END MODULE lmdz_fft_fftw
Note: See TracBrowser for help on using the repository browser.