source: LMDZ4/trunk/libf/filtrez/mod_filtre_fft.F90 @ 1068

Last change on this file since 1068 was 994, checked in by Laurent Fairhead, 16 years ago

Renommage filtre_fft.F90 en mod_filtre_fft.F90
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.7 KB
RevLine 
[994]1MODULE mod_filtre_fft
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,nlat,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) :: nlat
114    INTEGER,INTENT(IN) :: jj_begin
115    INTEGER,INTENT(IN) :: jj_end
116    INTEGER,INTENT(IN) :: nbniv
117    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
118
119    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
120    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
121    INTEGER            :: nb_vect
122    INTEGER :: i,j,l
123    INTEGER :: ll_nb
124!    REAL               :: vect_tmp(iim+inc,jj_end-jj_begin+1,nbniv)
125   
126    ll_nb=0
127!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
128    DO l=1,nbniv
129      ll_nb=ll_nb+1
130      DO j=1,jj_end-jj_begin+1
131        DO i=1,iim+1
132          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
133        ENDDO
134      ENDDO
135    ENDDO
136!$OMP END DO NOWAIT
137
138    nb_vect=(jj_end-jj_begin+1)*ll_nb
139
140!    vect_tmp=vect
141
142    CALL FFT_forward(vect,TF_vect,nb_vect)
143
144    DO l=1,ll_nb
145      DO j=1,jj_end-jj_begin+1
146        DO i=1,iim/2+1
147          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_u(i,jj_begin+j-1)
148        ENDDO
149      ENDDO
150    ENDDO
151 
152    CALL FFT_backward(TF_vect,vect,nb_vect)
153     
154!      PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx"
155!      DO j=1,jj_end-jj_begin+1
156!          DO i=1,iim
157!          PRINT *,"====",j,"----->",(vect_tmp(i,j,1)-vect(i,j,1))/ &
158!                                    ((vect_tmp(i,j,1)+vect(i,j,1))*0.5+1e-30)
159!          ENDDO
160!      ENDDO
161     
162    ll_nb=0
163!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
164    DO l=1,nbniv
165      ll_nb=ll_nb+1
166      DO j=1,jj_end-jj_begin+1
167        DO i=1,iim+1
168          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
169        ENDDO
170      ENDDO
171    ENDDO
172!$OMP END DO NOWAIT
173
174  END SUBROUTINE Filtre_u_fft
175 
176
177  SUBROUTINE Filtre_v_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
178    USE mod_fft
179#ifdef CPP_PARA
180    USE parallel,ONLY : OMP_CHUNK
181#endif
182    IMPLICIT NONE
183    INCLUDE 'dimensions.h'
184    INTEGER,INTENT(IN) :: nlat
185    INTEGER,INTENT(IN) :: jj_begin
186    INTEGER,INTENT(IN) :: jj_end
187    INTEGER,INTENT(IN) :: nbniv
188    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
189
190    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
191    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
192    INTEGER            :: nb_vect
193    INTEGER :: i,j,l
194    INTEGER :: ll_nb
195   
196    ll_nb=0
197!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
198    DO l=1,nbniv
199      ll_nb=ll_nb+1
200      DO j=1,jj_end-jj_begin+1
201        DO i=1,iim+1
202          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
203        ENDDO
204      ENDDO
205    ENDDO
206!$OMP END DO NOWAIT
207
208   
209    nb_vect=(jj_end-jj_begin+1)*ll_nb
210
211    CALL FFT_forward(vect,TF_vect,nb_vect)
212 
213    DO l=1,ll_nb
214      DO j=1,jj_end-jj_begin+1
215        DO i=1,iim/2+1
216          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_v(i,jj_begin+j-1)
217        ENDDO
218      ENDDO
219    ENDDO
220 
221    CALL FFT_backward(TF_vect,vect,nb_vect)
222   
223   
224    ll_nb=0
225!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
226    DO l=1,nbniv
227      ll_nb=ll_nb+1
228      DO j=1,jj_end-jj_begin+1
229        DO i=1,iim+1
230          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
231        ENDDO
232      ENDDO
233    ENDDO
234!$OMP END DO NOWAIT
235 
236  END SUBROUTINE Filtre_v_fft
237
238
239  SUBROUTINE Filtre_inv_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
240    USE mod_fft
241#ifdef CPP_PARA
242    USE parallel,ONLY : OMP_CHUNK
243#endif
244    IMPLICIT NONE
245    INCLUDE 'dimensions.h'
246    INTEGER,INTENT(IN) :: nlat
247    INTEGER,INTENT(IN) :: jj_begin
248    INTEGER,INTENT(IN) :: jj_end
249    INTEGER,INTENT(IN) :: nbniv
250    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
251
252     REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
253    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
254    INTEGER            :: nb_vect
255    INTEGER :: i,j,l
256    INTEGER :: ll_nb
257   
258    ll_nb=0
259!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
260    DO l=1,nbniv
261      ll_nb=ll_nb+1
262      DO j=1,jj_end-jj_begin+1
263        DO i=1,iim+1
264          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
265        ENDDO
266      ENDDO
267    ENDDO
268!$OMP END DO NOWAIT
269
270    nb_vect=(jj_end-jj_begin+1)*ll_nb
271
272    CALL FFT_forward(vect,TF_vect,nb_vect)
273 
274    DO l=1,ll_nb
275      DO j=1,jj_end-jj_begin+1
276        DO i=1,iim/2+1
277          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_inv(i,jj_begin+j-1)
278        ENDDO
279      ENDDO
280    ENDDO
281 
282    CALL FFT_backward(TF_vect,vect,nb_vect)
283
284    ll_nb=0
285!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
286    DO l=1,nbniv
287      ll_nb=ll_nb+1
288      DO j=1,jj_end-jj_begin+1
289        DO i=1,iim+1
290          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
291        ENDDO
292      ENDDO
293    ENDDO
294!$OMP END DO NOWAIT
295
296  END SUBROUTINE Filtre_inv_fft 
297 
298 
299!  SUBROUTINE get_ll_index(nbniv,ll_index,ll_nb)
300!  IMPLICIT NONE
301!    INTEGER,INTENT(IN)  :: nbniv
302!    INTEGER,INTENT(OUT) :: ll_index(nbniv)
303!    INTEGER,INTENT(OUT) :: ll_nb
304!
305!    INTEGER :: l,ll_begin, ll_end
306!   INTEGER :: omp_rank,omp_size
307!   INTEGER :: OMP_GET_NUM_THREADS
308!   INTEGER :: omp_chunk
309!   EXTERNAL OMP_GET_NUM_THREADS
310!   INTEGER :: OMP_GET_THREAD_NUM
311!   EXTERNAL OMP_GET_THREAD_NUM
312!
313!   
314!   omp_size=OMP_GET_NUM_THREADS()
315!   omp_rank=OMP_GET_THREAD_NUM()   
316!   omp_chunk=nbniv/omp_size+min(1,MOD(nbniv,omp_size))
317!   
318!   ll_begin=omp_rank*OMP_CHUNK+1
319!   ll_nb=0
320!   DO WHILE (ll_begin<=nbniv)
321!     ll_end=min(ll_begin+OMP_CHUNK-1,nbniv)
322!     DO l=ll_begin,ll_end
323!       ll_nb=ll_nb+1
324!       ll_index(ll_nb)=l
325!     ENDDO
326!     ll_begin=ll_begin+omp_size*OMP_CHUNK
327!   ENDDO
328! 
329!  END SUBROUTINE get_ll_index
330   
331END MODULE mod_filtre_fft
332 
Note: See TracBrowser for help on using the repository browser.