My Project
 All Classes Files Functions Variables Macros
mod_filtre_fft.F90
Go to the documentation of this file.
1 !
2 ! $Id: mod_filtre_fft.F90 1403 2010-07-01 09:02:53Z fairhead $
3 !
4 
5 MODULE mod_filtre_fft
6 
7  LOGICAL,SAVE :: use_filtre_fft
8  REAL,SAVE,ALLOCATABLE :: Filtre_u(:,:)
9  REAL,SAVE,ALLOCATABLE :: Filtre_v(:,:)
10  REAL,SAVE,ALLOCATABLE :: Filtre_inv(:,:)
11 
12 CONTAINS
13 
14  SUBROUTINE init_filtre_fft(coeffu,modfrstu,jfiltnu,jfiltsu,coeffv,modfrstv,jfiltnv,jfiltsv)
15  USE mod_fft
16  IMPLICIT NONE
17  include 'dimensions.h'
18  REAL, INTENT(IN) :: coeffu(iim,jjm)
19  INTEGER,INTENT(IN) :: modfrstu(jjm)
20  INTEGER,INTENT(IN) :: jfiltnu
21  INTEGER,INTENT(IN) :: jfiltsu
22  REAL, INTENT(IN) :: coeffv(iim,jjm)
23  INTEGER,INTENT(IN) :: modfrstv(jjm)
24  INTEGER,INTENT(IN) :: jfiltnv
25  INTEGER,INTENT(IN) :: jfiltsv
26 
27  INTEGER :: index_vp(iim)
28  INTEGER :: i,j
29  INTEGER :: l,ll_nb
30 
31  index_vp(1)=1
32  DO i=1,iim/2
33  index_vp(i+1)=i*2
34  ENDDO
35 
36  DO i=1,iim/2-1
37  index_vp(iim/2+i+1)=iim-2*i+1
38  ENDDO
39 
40  ALLOCATE(filtre_u(iim,jjm))
41  ALLOCATE(filtre_v(iim,jjm))
42  ALLOCATE(filtre_inv(iim,jjm))
43 
44 
45  DO j=2,jfiltnu
46  DO i=1,iim
47  IF (index_vp(i) < modfrstu(j)) THEN
48  filtre_u(i,j)=0
49  ELSE
50  filtre_u(i,j)=coeffu(index_vp(i),j)
51  ENDIF
52  ENDDO
53  ENDDO
54 
55  DO j=jfiltsu,jjm
56  DO i=1,iim
57  IF (index_vp(i) < modfrstu(j)) THEN
58  filtre_u(i,j)=0
59  ELSE
60  filtre_u(i,j)=coeffu(index_vp(i),j)
61  ENDIF
62  ENDDO
63  ENDDO
64 
65  DO j=1,jfiltnv
66  DO i=1,iim
67  IF (index_vp(i) < modfrstv(j)) THEN
68  filtre_v(i,j)=0
69  ELSE
70  filtre_v(i,j)=coeffv(index_vp(i),j)
71  ENDIF
72  ENDDO
73  ENDDO
74 
75  DO j=jfiltsv,jjm
76  DO i=1,iim
77  IF (index_vp(i) < modfrstv(j)) THEN
78  filtre_v(i,j)=0
79  ELSE
80  filtre_v(i,j)=coeffv(index_vp(i),j)
81  ENDIF
82  ENDDO
83  ENDDO
84 
85  DO j=2,jfiltnu
86  DO i=1,iim
87  IF (index_vp(i) < modfrstu(j)) THEN
88  filtre_inv(i,j)=0
89  ELSE
90  filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
91  ENDIF
92  ENDDO
93  ENDDO
94 
95  DO j=jfiltsu,jjm
96  DO i=1,iim
97  IF (index_vp(i) < modfrstu(j)) THEN
98  filtre_inv(i,j)=0
99  ELSE
100  filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
101  ENDIF
102  ENDDO
103  ENDDO
104 
105 #ifdef FFT_FFTW
106 
107  WRITE (*,*)"COTH jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv"
108  WRITE (*,*)jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv
109  WRITE (*,*)max(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1
110  CALL init_fft(iim,(llm+1)*(max(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1))
111 #else
112  CALL init_fft(iim,(jjm+1)*(llm+1))
113 #endif
114 
115  END SUBROUTINE init_filtre_fft
116 
117  SUBROUTINE filtre_u_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
118  USE mod_fft
119 #ifdef CPP_PARA
120  USE parallel,ONLY : OMP_CHUNK
121 #endif
122  IMPLICIT NONE
123  include 'dimensions.h'
124  INTEGER,INTENT(IN) :: nlat
125  INTEGER,INTENT(IN) :: jj_begin
126  INTEGER,INTENT(IN) :: jj_end
127  INTEGER,INTENT(IN) :: nbniv
128  REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
129 
130  REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
131  COMPLEX :: tf_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
132  INTEGER :: nb_vect
133  INTEGER :: i,j,l
134  INTEGER :: ll_nb
135 
136  ll_nb=0
137 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
138  DO l=1,nbniv
139  ll_nb=ll_nb+1
140  DO j=1,jj_end-jj_begin+1
141  DO i=1,iim+1
142  vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
143  ENDDO
144  ENDDO
145  ENDDO
146 !$OMP END DO NOWAIT
147 
148  nb_vect=(jj_end-jj_begin+1)*ll_nb
149 
150  CALL fft_forward(vect,tf_vect,nb_vect)
151 
152  DO l=1,ll_nb
153  DO j=1,jj_end-jj_begin+1
154  DO i=1,iim/2+1
155  tf_vect(i,j,l)=tf_vect(i,j,l)*filtre_u(i,jj_begin+j-1)
156  ENDDO
157  ENDDO
158  ENDDO
159 
160  CALL fft_backward(tf_vect,vect,nb_vect)
161 
162 
163  ll_nb=0
164 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
165  DO l=1,nbniv
166  ll_nb=ll_nb+1
167  DO j=1,jj_end-jj_begin+1
168  DO i=1,iim+1
169  vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
170  ENDDO
171  ENDDO
172  ENDDO
173 !$OMP END DO NOWAIT
174 
175  END SUBROUTINE filtre_u_fft
176 
177 
178  SUBROUTINE filtre_v_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
179  USE mod_fft
180 #ifdef CPP_PARA
181  USE parallel,ONLY : OMP_CHUNK
182 #endif
183  IMPLICIT NONE
184  include 'dimensions.h'
185  INTEGER,INTENT(IN) :: nlat
186  INTEGER,INTENT(IN) :: jj_begin
187  INTEGER,INTENT(IN) :: jj_end
188  INTEGER,INTENT(IN) :: nbniv
189  REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
190 
191  REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
192  COMPLEX :: tf_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
193  INTEGER :: nb_vect
194  INTEGER :: i,j,l
195  INTEGER :: ll_nb
196 
197  ll_nb=0
198 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
199  DO l=1,nbniv
200  ll_nb=ll_nb+1
201  DO j=1,jj_end-jj_begin+1
202  DO i=1,iim+1
203  vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
204  ENDDO
205  ENDDO
206  ENDDO
207 !$OMP END DO NOWAIT
208 
209 
210  nb_vect=(jj_end-jj_begin+1)*ll_nb
211 
212  CALL fft_forward(vect,tf_vect,nb_vect)
213 
214  DO l=1,ll_nb
215  DO j=1,jj_end-jj_begin+1
216  DO i=1,iim/2+1
217  tf_vect(i,j,l)=tf_vect(i,j,l)*filtre_v(i,jj_begin+j-1)
218  ENDDO
219  ENDDO
220  ENDDO
221 
222  CALL fft_backward(tf_vect,vect,nb_vect)
223 
224 
225  ll_nb=0
226 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
227  DO l=1,nbniv
228  ll_nb=ll_nb+1
229  DO j=1,jj_end-jj_begin+1
230  DO i=1,iim+1
231  vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
232  ENDDO
233  ENDDO
234  ENDDO
235 !$OMP END DO NOWAIT
236 
237  END SUBROUTINE filtre_v_fft
238 
239 
240  SUBROUTINE filtre_inv_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
241  USE mod_fft
242 #ifdef CPP_PARA
243  USE parallel,ONLY : OMP_CHUNK
244 #endif
245  IMPLICIT NONE
246  include 'dimensions.h'
247  INTEGER,INTENT(IN) :: nlat
248  INTEGER,INTENT(IN) :: jj_begin
249  INTEGER,INTENT(IN) :: jj_end
250  INTEGER,INTENT(IN) :: nbniv
251  REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
252 
253  REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
254  COMPLEX :: tf_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
255  INTEGER :: nb_vect
256  INTEGER :: i,j,l
257  INTEGER :: ll_nb
258 
259  ll_nb=0
260 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
261  DO l=1,nbniv
262  ll_nb=ll_nb+1
263  DO j=1,jj_end-jj_begin+1
264  DO i=1,iim+1
265  vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
266  ENDDO
267  ENDDO
268  ENDDO
269 !$OMP END DO NOWAIT
270 
271  nb_vect=(jj_end-jj_begin+1)*ll_nb
272 
273  CALL fft_forward(vect,tf_vect,nb_vect)
274 
275  DO l=1,ll_nb
276  DO j=1,jj_end-jj_begin+1
277  DO i=1,iim/2+1
278  tf_vect(i,j,l)=tf_vect(i,j,l)*filtre_inv(i,jj_begin+j-1)
279  ENDDO
280  ENDDO
281  ENDDO
282 
283  CALL fft_backward(tf_vect,vect,nb_vect)
284 
285  ll_nb=0
286 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
287  DO l=1,nbniv
288  ll_nb=ll_nb+1
289  DO j=1,jj_end-jj_begin+1
290  DO i=1,iim+1
291  vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
292  ENDDO
293  ENDDO
294  ENDDO
295 !$OMP END DO NOWAIT
296 
297  END SUBROUTINE filtre_inv_fft
298 
299 END MODULE mod_filtre_fft
300