source: LMDZ5/trunk/libf/filtrez/mod_filtre_fft_loc.F90 @ 1679

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

Module manquant pour dynamique localisée

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