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

Last change on this file since 1373 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.1 KB
Line 
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!    REAL               :: vect_test(iim+inc,jj_end-jj_begin+1,nbniv)
121    COMPLEX*16         :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
122!    COMPLEX*16         :: TF_vect_test(iim/2+1,jj_end-jj_begin+1,nbniv)
123    INTEGER            :: nb_vect
124    INTEGER :: i,j,l
125    INTEGER :: ll_nb
126!    REAL               :: vect_tmp(iim+inc,jj_end-jj_begin+1,nbniv)
127   
128    ll_nb=0
129!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
130    DO l=1,nbniv
131      ll_nb=ll_nb+1
132      DO j=1,jj_end-jj_begin+1
133        DO i=1,iim+1
134          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
135        ENDDO
136      ENDDO
137    ENDDO
138!$OMP END DO NOWAIT
139
140    nb_vect=(jj_end-jj_begin+1)*ll_nb
141
142!    vect_tmp=vect
143
144    CALL FFT_forward(vect,TF_vect,nb_vect)
145
146!    CALL FFT_forward(vect,TF_vect_test,nb_vect)
147!      PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx"
148!      DO j=1,jj_end-jj_begin+1
149!      DO i=1,iim/2+1
150!         PRINT *,"====",i,j,"----->",TF_vect_test(i,j,1)
151!       ENDDO
152!      ENDDO
153
154    DO l=1,ll_nb
155      DO j=1,jj_end-jj_begin+1
156        DO i=1,iim/2+1
157          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_u(i,jj_begin+j-1)
158        ENDDO
159      ENDDO
160    ENDDO
161       
162    CALL FFT_backward(TF_vect,vect,nb_vect)
163!    CALL FFT_backward(TF_vect_test,vect_test,nb_vect)
164         
165!      PRINT *,"XXXXXXXXXXXXX Filtre_u_FFT xxxxxxxxxxxx"
166!      DO j=1,jj_end-jj_begin+1
167!         DO i=1,iim
168!           PRINT *,"====",i,j,"----->",vect_test(i,j,1)
169!         ENDDO
170!      ENDDO
171
172    ll_nb=0
173!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
174    DO l=1,nbniv
175      ll_nb=ll_nb+1
176      DO j=1,jj_end-jj_begin+1
177        DO i=1,iim+1
178          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
179        ENDDO
180      ENDDO
181    ENDDO
182!$OMP END DO NOWAIT
183
184  END SUBROUTINE Filtre_u_fft
185 
186
187  SUBROUTINE Filtre_v_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
188    USE mod_fft
189#ifdef CPP_PARA
190    USE parallel,ONLY : OMP_CHUNK
191#endif
192    IMPLICIT NONE
193    INCLUDE 'dimensions.h'
194    INTEGER,INTENT(IN) :: nlat
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,nlat,nbniv)
199
200    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
201    COMPLEX*16         :: 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,nlat,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) :: nlat
257    INTEGER,INTENT(IN) :: jj_begin
258    INTEGER,INTENT(IN) :: jj_end
259    INTEGER,INTENT(IN) :: nbniv
260    REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
261
262    REAL               :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
263    COMPLEX*16         :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
264    INTEGER            :: nb_vect
265    INTEGER :: i,j,l
266    INTEGER :: ll_nb
267   
268    ll_nb=0
269!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
270    DO l=1,nbniv
271      ll_nb=ll_nb+1
272      DO j=1,jj_end-jj_begin+1
273        DO i=1,iim+1
274          vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
275        ENDDO
276      ENDDO
277    ENDDO
278!$OMP END DO NOWAIT
279
280    nb_vect=(jj_end-jj_begin+1)*ll_nb
281
282    CALL FFT_forward(vect,TF_vect,nb_vect)
283 
284    DO l=1,ll_nb
285      DO j=1,jj_end-jj_begin+1
286        DO i=1,iim/2+1
287          TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_inv(i,jj_begin+j-1)
288        ENDDO
289      ENDDO
290    ENDDO
291 
292    CALL FFT_backward(TF_vect,vect,nb_vect)
293
294    ll_nb=0
295!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
296    DO l=1,nbniv
297      ll_nb=ll_nb+1
298      DO j=1,jj_end-jj_begin+1
299        DO i=1,iim+1
300          vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
301        ENDDO
302      ENDDO
303    ENDDO
304!$OMP END DO NOWAIT
305
306  END SUBROUTINE Filtre_inv_fft 
307 
308 
309!  SUBROUTINE get_ll_index(nbniv,ll_index,ll_nb)
310!  IMPLICIT NONE
311!    INTEGER,INTENT(IN)  :: nbniv
312!    INTEGER,INTENT(OUT) :: ll_index(nbniv)
313!    INTEGER,INTENT(OUT) :: ll_nb
314!
315!    INTEGER :: l,ll_begin, ll_end
316!   INTEGER :: omp_rank,omp_size
317!   INTEGER :: OMP_GET_NUM_THREADS
318!   INTEGER :: omp_chunk
319!   EXTERNAL OMP_GET_NUM_THREADS
320!   INTEGER :: OMP_GET_THREAD_NUM
321!   EXTERNAL OMP_GET_THREAD_NUM
322!
323!   
324!   omp_size=OMP_GET_NUM_THREADS()
325!   omp_rank=OMP_GET_THREAD_NUM()   
326!   omp_chunk=nbniv/omp_size+min(1,MOD(nbniv,omp_size))
327!   
328!   ll_begin=omp_rank*OMP_CHUNK+1
329!   ll_nb=0
330!   DO WHILE (ll_begin<=nbniv)
331!     ll_end=min(ll_begin+OMP_CHUNK-1,nbniv)
332!     DO l=ll_begin,ll_end
333!       ll_nb=ll_nb+1
334!       ll_index(ll_nb)=l
335!     ENDDO
336!     ll_begin=ll_begin+omp_size*OMP_CHUNK
337!   ENDDO
338! 
339!  END SUBROUTINE get_ll_index
340   
341END MODULE mod_filtre_fft
342 
Note: See TracBrowser for help on using the repository browser.