Context Navigation

mod_filtre_fft_loc.F90 @ 5308

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

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_lmdz,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_lmdz,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_lmdz,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