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

Last change on this file since 3906 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.
  • 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
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
[1823]109    USE parallel_lmdz,ONLY : OMP_CHUNK
[1679]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)
[1701]121    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
[1679]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
[1823]189    USE parallel_lmdz,ONLY : OMP_CHUNK
[1679]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)
[1701]201    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
[1679]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
[1823]252    USE parallel_lmdz,ONLY : OMP_CHUNK
[1679]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)
[1701]264    COMPLEX            :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
[1679]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.