source: LMDZ6/trunk/libf/filtrez/mod_filtre_fft_loc.F90 @ 5295

Last change on this file since 5295 was 5271, checked in by abarral, 8 days ago

Move dimensions.h into a module
Nb: doesn't compile yet
  • 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
File size: 8.1 KB
RevLine 
[1679]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
[5271]12    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
13IMPLICIT NONE
14
[1679]15    REAL,   INTENT(IN) :: coeffu(iim,jjm)
16    INTEGER,INTENT(IN) :: modfrstu(jjm)
17    INTEGER,INTENT(IN) :: jfiltnu
18    INTEGER,INTENT(IN) :: jfiltsu
19    REAL,   INTENT(IN) :: coeffv(iim,jjm)
20    INTEGER,INTENT(IN) :: modfrstv(jjm)
21    INTEGER,INTENT(IN) :: jfiltnv
22    INTEGER,INTENT(IN) :: jfiltsv
23   
24    INTEGER            :: index_vp(iim)
25    INTEGER            :: i,j
26   
27    index_vp(1)=1
28    DO i=1,iim/2
29      index_vp(i+1)=i*2
30    ENDDO
31   
32    DO i=1,iim/2-1
33      index_vp(iim/2+i+1)=iim-2*i+1
34    ENDDO
35   
36    ALLOCATE(Filtre_u(iim,jjm))
37    ALLOCATE(Filtre_v(iim,jjm))
38    ALLOCATE(Filtre_inv(iim,jjm))
39 
40   
41    DO j=2,jfiltnu
42      DO i=1,iim
43        IF (index_vp(i) < modfrstu(j)) THEN
44          Filtre_u(i,j)=0
45        ELSE
46          Filtre_u(i,j)=coeffu(index_vp(i),j)
47        ENDIF
48      ENDDO
49    ENDDO
50   
51    DO j=jfiltsu,jjm
52      DO i=1,iim
53        IF (index_vp(i) < modfrstu(j)) THEN
54          Filtre_u(i,j)=0
55        ELSE
56          Filtre_u(i,j)=coeffu(index_vp(i),j)
57        ENDIF
58      ENDDO
59    ENDDO
60 
61    DO j=1,jfiltnv
62      DO i=1,iim
63        IF (index_vp(i) < modfrstv(j)) THEN
64          Filtre_v(i,j)=0
65        ELSE
66          Filtre_v(i,j)=coeffv(index_vp(i),j)
67        ENDIF
68      ENDDO
69    ENDDO
70   
71    DO j=jfiltsv,jjm
72      DO i=1,iim
73        IF (index_vp(i) < modfrstv(j)) THEN
74          Filtre_v(i,j)=0
75        ELSE
76          Filtre_v(i,j)=coeffv(index_vp(i),j)
77        ENDIF
78      ENDDO
79    ENDDO
80         
81    DO j=2,jfiltnu
82      DO i=1,iim
83        IF (index_vp(i) < modfrstu(j)) THEN
84          Filtre_inv(i,j)=0
85        ELSE
86          Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
87        ENDIF
88      ENDDO
89    ENDDO
90
91    DO j=jfiltsu,jjm
92      DO i=1,iim
93        IF (index_vp(i) < modfrstu(j)) THEN
94          Filtre_inv(i,j)=0
95        ELSE
96          Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
97        ENDIF
98      ENDDO
99    ENDDO
100   
101   
102!    CALL Init_FFT(iim,(jjm+1)*(llm+1))
103       
104   
105  END SUBROUTINE Init_filtre_fft
106 
107  SUBROUTINE Filtre_u_fft(vect_inout,jjb,jje,jj_begin,jj_end,nbniv)
108    USE mod_fft
109#ifdef CPP_PARA
[1823]110    USE parallel_lmdz,ONLY : OMP_CHUNK
[1679]111#endif
[5271]112    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
113IMPLICIT NONE
114
[1679]115    INTEGER,INTENT(IN) :: jjb
116    INTEGER,INTENT(IN) :: jje
117    INTEGER,INTENT(IN) :: jj_begin
118    INTEGER,INTENT(IN) :: jj_end
119    INTEGER,INTENT(IN) :: nbniv
120    REAL,INTENT(INOUT) :: vect_inout(iim+1,jjb:jje,nbniv)
121
122    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
[1701]123    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
[1679]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
[1823]191    USE parallel_lmdz,ONLY : OMP_CHUNK
[1679]192#endif
[5271]193    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
194IMPLICIT NONE
195
[1679]196    INTEGER,INTENT(IN) :: jjb
197    INTEGER,INTENT(IN) :: jje
198    INTEGER,INTENT(IN) :: jj_begin
199    INTEGER,INTENT(IN) :: jj_end
200    INTEGER,INTENT(IN) :: nbniv
201    REAL,INTENT(INOUT) :: vect_inout(iim+1,jjb:jje,nbniv)
202
203    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
[1701]204    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
[1679]205    INTEGER            :: nb_vect
206    INTEGER :: i,j,l
207    INTEGER :: ll_nb
208   
209    ll_nb=0
210!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
211    DO l=1,nbniv
212      ll_nb=ll_nb+1
213      DO j=1,jj_end-jj_begin+1
214        DO i=1,iim+1
215          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
216        ENDDO
217      ENDDO
218    ENDDO
219!$OMP END DO NOWAIT
220
221   
222    nb_vect=(jj_end-jj_begin+1)*ll_nb
223
224    CALL FFT_forward(vect,TF_vect,nb_vect)
225 
226    DO l=1,ll_nb
227      DO j=1,jj_end-jj_begin+1
228        DO i=1,iim/2+1
229          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_v(i,jj_begin+j-1)
230        ENDDO
231      ENDDO
232    ENDDO
233 
234    CALL FFT_backward(TF_vect,vect,nb_vect)
235   
236   
237    ll_nb=0
238!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
239    DO l=1,nbniv
240      ll_nb=ll_nb+1
241      DO j=1,jj_end-jj_begin+1
242        DO i=1,iim+1
243          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
244        ENDDO
245      ENDDO
246    ENDDO
247!$OMP END DO NOWAIT
248 
249  END SUBROUTINE Filtre_v_fft
250
251
252  SUBROUTINE Filtre_inv_fft(vect_inout,jjb,jje,jj_begin,jj_end,nbniv)
253    USE mod_fft
254#ifdef CPP_PARA
[1823]255    USE parallel_lmdz,ONLY : OMP_CHUNK
[1679]256#endif
[5271]257    USE dimensions_mod, ONLY: iim, jjm, llm, ndm
258IMPLICIT NONE
259
[1679]260    INTEGER,INTENT(IN) :: jjb
261    INTEGER,INTENT(IN) :: jje
262    INTEGER,INTENT(IN) :: jj_begin
263    INTEGER,INTENT(IN) :: jj_end
264    INTEGER,INTENT(IN) :: nbniv
265    REAL,INTENT(INOUT) :: vect_inout(iim+1,jjb:jje,nbniv)
266
267    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
[1701]268    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
[1679]269    INTEGER            :: nb_vect
270    INTEGER :: i,j,l
271    INTEGER :: ll_nb
272   
273    ll_nb=0
274!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
275    DO l=1,nbniv
276      ll_nb=ll_nb+1
277      DO j=1,jj_end-jj_begin+1
278        DO i=1,iim+1
279          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
280        ENDDO
281      ENDDO
282    ENDDO
283!$OMP END DO NOWAIT
284
285    nb_vect=(jj_end-jj_begin+1)*ll_nb
286
287    CALL FFT_forward(vect,TF_vect,nb_vect)
288 
289    DO l=1,ll_nb
290      DO j=1,jj_end-jj_begin+1
291        DO i=1,iim/2+1
292          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_inv(i,jj_begin+j-1)
293        ENDDO
294      ENDDO
295    ENDDO
296 
297    CALL FFT_backward(TF_vect,vect,nb_vect)
298
299    ll_nb=0
300!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
301    DO l=1,nbniv
302      ll_nb=ll_nb+1
303      DO j=1,jj_end-jj_begin+1
304        DO i=1,iim+1
305          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
306        ENDDO
307      ENDDO
308    ENDDO
309!$OMP END DO NOWAIT
310
311  END SUBROUTINE Filtre_inv_fft 
312 
313 
314!  SUBROUTINE get_ll_index(nbniv,ll_index,ll_nb)
315!  IMPLICIT NONE
316!    INTEGER,INTENT(IN)  :: nbniv
317!    INTEGER,INTENT(OUT) :: ll_index(nbniv)
318!    INTEGER,INTENT(OUT) :: ll_nb
319!
320!    INTEGER :: l,ll_begin, ll_end
321!   INTEGER :: omp_rank,omp_size
322!   INTEGER :: OMP_GET_NUM_THREADS
323!   INTEGER :: omp_chunk
324!   EXTERNAL OMP_GET_NUM_THREADS
325!   INTEGER :: OMP_GET_THREAD_NUM
326!   EXTERNAL OMP_GET_THREAD_NUM
327!
328!   
329!   omp_size=OMP_GET_NUM_THREADS()
330!   omp_rank=OMP_GET_THREAD_NUM()   
331!   omp_chunk=nbniv/omp_size+min(1,MOD(nbniv,omp_size))
332!   
333!   ll_begin=omp_rank*OMP_CHUNK+1
334!   ll_nb=0
335!   DO WHILE (ll_begin<=nbniv)
336!     ll_end=min(ll_begin+OMP_CHUNK-1,nbniv)
337!     DO l=ll_begin,ll_end
338!       ll_nb=ll_nb+1
339!       ll_index(ll_nb)=l
340!     ENDDO
341!     ll_begin=ll_begin+omp_size*OMP_CHUNK
342!   ENDDO
343! 
344!  END SUBROUTINE get_ll_index
345   
346END MODULE mod_filtre_fft_loc
347 
Note: See TracBrowser for help on using the repository browser.