source: LMDZ5/branches/testing/libf/filtrez/mod_filtre_fft_loc.F90 @ 1707

Last change on this file since 1707 was 1707, checked in by Laurent Fairhead, 12 years ago

Version testing basée sur la r1706


Testing release based on r1706

File size: 8.1 KB
Line 
1MODULE mod_filtre_fft_loc
2
3  LOGICAL,SAVE :: use_filtre_fft
4  REAL,SAVE,ALLOCATABLE :: Filtre_u(:,:)
5  REAL,SAVE,ALLOCATABLE :: Filtre_v(:,:)
6  REAL,SAVE,ALLOCATABLE :: Filtre_inv(:,:)
7
8CONTAINS
9 
10  SUBROUTINE Init_filtre_fft(coeffu,modfrstu,jfiltnu,jfiltsu,coeffv,modfrstv,jfiltnv,jfiltsv)
11    USE mod_fft
12    IMPLICIT NONE
13    include 'dimensions.h'
14    REAL,   INTENT(IN) :: coeffu(iim,jjm)
15    INTEGER,INTENT(IN) :: modfrstu(jjm)
16    INTEGER,INTENT(IN) :: jfiltnu
17    INTEGER,INTENT(IN) :: jfiltsu
18    REAL,   INTENT(IN) :: coeffv(iim,jjm)
19    INTEGER,INTENT(IN) :: modfrstv(jjm)
20    INTEGER,INTENT(IN) :: jfiltnv
21    INTEGER,INTENT(IN) :: jfiltsv
22   
23    INTEGER            :: index_vp(iim)
24    INTEGER            :: i,j
25   
26    index_vp(1)=1
27    DO i=1,iim/2
28      index_vp(i+1)=i*2
29    ENDDO
30   
31    DO i=1,iim/2-1
32      index_vp(iim/2+i+1)=iim-2*i+1
33    ENDDO
34   
35    ALLOCATE(Filtre_u(iim,jjm))
36    ALLOCATE(Filtre_v(iim,jjm))
37    ALLOCATE(Filtre_inv(iim,jjm))
38 
39   
40    DO j=2,jfiltnu
41      DO i=1,iim
42        IF (index_vp(i) < modfrstu(j)) THEN
43          Filtre_u(i,j)=0
44        ELSE
45          Filtre_u(i,j)=coeffu(index_vp(i),j)
46        ENDIF
47      ENDDO
48    ENDDO
49   
50    DO j=jfiltsu,jjm
51      DO i=1,iim
52        IF (index_vp(i) < modfrstu(j)) THEN
53          Filtre_u(i,j)=0
54        ELSE
55          Filtre_u(i,j)=coeffu(index_vp(i),j)
56        ENDIF
57      ENDDO
58    ENDDO
59 
60    DO j=1,jfiltnv
61      DO i=1,iim
62        IF (index_vp(i) < modfrstv(j)) THEN
63          Filtre_v(i,j)=0
64        ELSE
65          Filtre_v(i,j)=coeffv(index_vp(i),j)
66        ENDIF
67      ENDDO
68    ENDDO
69   
70    DO j=jfiltsv,jjm
71      DO i=1,iim
72        IF (index_vp(i) < modfrstv(j)) THEN
73          Filtre_v(i,j)=0
74        ELSE
75          Filtre_v(i,j)=coeffv(index_vp(i),j)
76        ENDIF
77      ENDDO
78    ENDDO
79         
80    DO j=2,jfiltnu
81      DO i=1,iim
82        IF (index_vp(i) < modfrstu(j)) THEN
83          Filtre_inv(i,j)=0
84        ELSE
85          Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
86        ENDIF
87      ENDDO
88    ENDDO
89
90    DO j=jfiltsu,jjm
91      DO i=1,iim
92        IF (index_vp(i) < modfrstu(j)) THEN
93          Filtre_inv(i,j)=0
94        ELSE
95          Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
96        ENDIF
97      ENDDO
98    ENDDO
99   
100   
101!    CALL Init_FFT(iim,(jjm+1)*(llm+1))
102       
103   
104  END SUBROUTINE Init_filtre_fft
105 
106  SUBROUTINE Filtre_u_fft(vect_inout,jjb,jje,jj_begin,jj_end,nbniv)
107    USE mod_fft
108#ifdef CPP_PARA
109    USE parallel,ONLY : OMP_CHUNK
110#endif
111    IMPLICIT NONE
112    include 'dimensions.h'
113    INTEGER,INTENT(IN) :: jjb
114    INTEGER,INTENT(IN) :: jje
115    INTEGER,INTENT(IN) :: jj_begin
116    INTEGER,INTENT(IN) :: jj_end
117    INTEGER,INTENT(IN) :: nbniv
118    REAL,INTENT(INOUT) :: vect_inout(iim+1,jjb:jje,nbniv)
119
120    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
121    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
122    INTEGER            :: nb_vect
123    INTEGER :: i,j,l
124    INTEGER :: ll_nb
125!    REAL               :: vect_tmp(iim+inc,jj_end-jj_begin+1,nbniv)
126   
127    ll_nb=0
128!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
129    DO l=1,nbniv
130      ll_nb=ll_nb+1
131      DO j=1,jj_end-jj_begin+1
132        DO i=1,iim+1
133          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
134        ENDDO
135      ENDDO
136    ENDDO
137!$OMP END DO NOWAIT
138
139    nb_vect=(jj_end-jj_begin+1)*ll_nb
140
141!    vect_tmp=vect
142
143    CALL FFT_forward(vect,TF_vect,nb_vect)
144
145!    CALL FFT_forward(vect,TF_vect_test,nb_vect)
146!      PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx"
147!      DO j=1,jj_end-jj_begin+1
148!      DO i=1,iim/2+1
149!         PRINT *,"====",i,j,"----->",TF_vect_test(i,j,1)
150!       ENDDO
151!      ENDDO
152
153    DO l=1,ll_nb
154      DO j=1,jj_end-jj_begin+1
155        DO i=1,iim/2+1
156          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_u(i,jj_begin+j-1)
157        ENDDO
158      ENDDO
159    ENDDO
160       
161    CALL FFT_backward(TF_vect,vect,nb_vect)
162!    CALL FFT_backward(TF_vect_test,vect_test,nb_vect)
163         
164!      PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx"
165!      DO j=1,jj_end-jj_begin+1
166!         DO i=1,iim
167!           PRINT *,"====",i,j,"----->",vect_test(i,j,1)
168!         ENDDO
169!      ENDDO
170
171    ll_nb=0
172!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
173    DO l=1,nbniv
174      ll_nb=ll_nb+1
175      DO j=1,jj_end-jj_begin+1
176        DO i=1,iim+1
177          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
178        ENDDO
179      ENDDO
180    ENDDO
181!$OMP END DO NOWAIT
182
183  END SUBROUTINE Filtre_u_fft
184 
185
186  SUBROUTINE Filtre_v_fft(vect_inout,jjb,jje,jj_begin,jj_end,nbniv)
187    USE mod_fft
188#ifdef CPP_PARA
189    USE parallel,ONLY : OMP_CHUNK
190#endif
191    IMPLICIT NONE
192    INCLUDE 'dimensions.h'
193    INTEGER,INTENT(IN) :: jjb
194    INTEGER,INTENT(IN) :: jje
195    INTEGER,INTENT(IN) :: jj_begin
196    INTEGER,INTENT(IN) :: jj_end
197    INTEGER,INTENT(IN) :: nbniv
198    REAL,INTENT(INOUT) :: vect_inout(iim+1,jjb:jje,nbniv)
199
200    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
201    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
202    INTEGER            :: nb_vect
203    INTEGER :: i,j,l
204    INTEGER :: ll_nb
205   
206    ll_nb=0
207!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
208    DO l=1,nbniv
209      ll_nb=ll_nb+1
210      DO j=1,jj_end-jj_begin+1
211        DO i=1,iim+1
212          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
213        ENDDO
214      ENDDO
215    ENDDO
216!$OMP END DO NOWAIT
217
218   
219    nb_vect=(jj_end-jj_begin+1)*ll_nb
220
221    CALL FFT_forward(vect,TF_vect,nb_vect)
222 
223    DO l=1,ll_nb
224      DO j=1,jj_end-jj_begin+1
225        DO i=1,iim/2+1
226          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_v(i,jj_begin+j-1)
227        ENDDO
228      ENDDO
229    ENDDO
230 
231    CALL FFT_backward(TF_vect,vect,nb_vect)
232   
233   
234    ll_nb=0
235!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
236    DO l=1,nbniv
237      ll_nb=ll_nb+1
238      DO j=1,jj_end-jj_begin+1
239        DO i=1,iim+1
240          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
241        ENDDO
242      ENDDO
243    ENDDO
244!$OMP END DO NOWAIT
245 
246  END SUBROUTINE Filtre_v_fft
247
248
249  SUBROUTINE Filtre_inv_fft(vect_inout,jjb,jje,jj_begin,jj_end,nbniv)
250    USE mod_fft
251#ifdef CPP_PARA
252    USE parallel,ONLY : OMP_CHUNK
253#endif
254    IMPLICIT NONE
255    INCLUDE 'dimensions.h'
256    INTEGER,INTENT(IN) :: jjb
257    INTEGER,INTENT(IN) :: jje
258    INTEGER,INTENT(IN) :: jj_begin
259    INTEGER,INTENT(IN) :: jj_end
260    INTEGER,INTENT(IN) :: nbniv
261    REAL,INTENT(INOUT) :: vect_inout(iim+1,jjb:jje,nbniv)
262
263    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
264    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
265    INTEGER            :: nb_vect
266    INTEGER :: i,j,l
267    INTEGER :: ll_nb
268   
269    ll_nb=0
270!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
271    DO l=1,nbniv
272      ll_nb=ll_nb+1
273      DO j=1,jj_end-jj_begin+1
274        DO i=1,iim+1
275          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
276        ENDDO
277      ENDDO
278    ENDDO
279!$OMP END DO NOWAIT
280
281    nb_vect=(jj_end-jj_begin+1)*ll_nb
282
283    CALL FFT_forward(vect,TF_vect,nb_vect)
284 
285    DO l=1,ll_nb
286      DO j=1,jj_end-jj_begin+1
287        DO i=1,iim/2+1
288          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_inv(i,jj_begin+j-1)
289        ENDDO
290      ENDDO
291    ENDDO
292 
293    CALL FFT_backward(TF_vect,vect,nb_vect)
294
295    ll_nb=0
296!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
297    DO l=1,nbniv
298      ll_nb=ll_nb+1
299      DO j=1,jj_end-jj_begin+1
300        DO i=1,iim+1
301          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
302        ENDDO
303      ENDDO
304    ENDDO
305!$OMP END DO NOWAIT
306
307  END SUBROUTINE Filtre_inv_fft 
308 
309 
310!  SUBROUTINE get_ll_index(nbniv,ll_index,ll_nb)
311!  IMPLICIT NONE
312!    INTEGER,INTENT(IN)  :: nbniv
313!    INTEGER,INTENT(OUT) :: ll_index(nbniv)
314!    INTEGER,INTENT(OUT) :: ll_nb
315!
316!    INTEGER :: l,ll_begin, ll_end
317!   INTEGER :: omp_rank,omp_size
318!   INTEGER :: OMP_GET_NUM_THREADS
319!   INTEGER :: omp_chunk
320!   EXTERNAL OMP_GET_NUM_THREADS
321!   INTEGER :: OMP_GET_THREAD_NUM
322!   EXTERNAL OMP_GET_THREAD_NUM
323!
324!   
325!   omp_size=OMP_GET_NUM_THREADS()
326!   omp_rank=OMP_GET_THREAD_NUM()   
327!   omp_chunk=nbniv/omp_size+min(1,MOD(nbniv,omp_size))
328!   
329!   ll_begin=omp_rank*OMP_CHUNK+1
330!   ll_nb=0
331!   DO WHILE (ll_begin<=nbniv)
332!     ll_end=min(ll_begin+OMP_CHUNK-1,nbniv)
333!     DO l=ll_begin,ll_end
334!       ll_nb=ll_nb+1
335!       ll_index(ll_nb)=l
336!     ENDDO
337!     ll_begin=ll_begin+omp_size*OMP_CHUNK
338!   ENDDO
339
340!  END SUBROUTINE get_ll_index
341   
342END MODULE mod_filtre_fft_loc
343 
Note: See TracBrowser for help on using the repository browser.