Context Navigation

mod_filtre_fft_loc.F90 @ 1790

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

Version testing basée sur la r1706

Testing release based on r1706

File size: 8.1 KB

Line
1	MODULE 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
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,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	COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
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
189	USE parallel,ONLY : OMP_CHUNK
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)
201	COMPLEX :: 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,jjb,jje,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) :: 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)
264	COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
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
342	END MODULE mod_filtre_fft_loc
343

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

Download in other formats:

Original Format