Context Navigation

mod_filtre_fft.F90 @ 994

Last change on this file since 994 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

Line
1	MODULE 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
8	CONTAINS
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
331	END MODULE mod_filtre_fft
332

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format