Context Navigation

mod_filtre_fft.F90 @ 3325

Last change on this file since 3325 was 1910, checked in by Laurent Fairhead, 11 years ago
Merged trunk changes r1860:1909 into testing branch
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 Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 6.9 KB

Line
1	!
2	! $Id: mod_filtre_fft.F90 1910 2013-11-29 08:40:25Z musat $
3	!
4
5	MODULE mod_filtre_fft
6
7	LOGICAL,SAVE :: use_filtre_fft
8	REAL,SAVE,ALLOCATABLE :: Filtre_u(:,:)
9	REAL,SAVE,ALLOCATABLE :: Filtre_v(:,:)
10	REAL,SAVE,ALLOCATABLE :: Filtre_inv(:,:)
11
12	CONTAINS
13
14	SUBROUTINE Init_filtre_fft(coeffu,modfrstu,jfiltnu,jfiltsu,coeffv,modfrstv,jfiltnv,jfiltsv)
15	USE mod_fft
16	IMPLICIT NONE
17	include 'dimensions.h'
18	REAL, INTENT(IN) :: coeffu(iim,jjm)
19	INTEGER,INTENT(IN) :: modfrstu(jjm)
20	INTEGER,INTENT(IN) :: jfiltnu
21	INTEGER,INTENT(IN) :: jfiltsu
22	REAL, INTENT(IN) :: coeffv(iim,jjm)
23	INTEGER,INTENT(IN) :: modfrstv(jjm)
24	INTEGER,INTENT(IN) :: jfiltnv
25	INTEGER,INTENT(IN) :: jfiltsv
26
27	INTEGER :: index_vp(iim)
28	INTEGER :: i,j
29	INTEGER :: l,ll_nb
30
31	index_vp(1)=1
32	DO i=1,iim/2
33	index_vp(i+1)=i*2
34	ENDDO
35
36	DO i=1,iim/2-1
37	index_vp(iim/2+i+1)=iim-2*i+1
38	ENDDO
39
40	ALLOCATE(Filtre_u(iim,jjm))
41	ALLOCATE(Filtre_v(iim,jjm))
42	ALLOCATE(Filtre_inv(iim,jjm))
43
44
45	DO j=2,jfiltnu
46	DO i=1,iim
47	IF (index_vp(i) < modfrstu(j)) THEN
48	Filtre_u(i,j)=0
49	ELSE
50	Filtre_u(i,j)=coeffu(index_vp(i),j)
51	ENDIF
52	ENDDO
53	ENDDO
54
55	DO j=jfiltsu,jjm
56	DO i=1,iim
57	IF (index_vp(i) < modfrstu(j)) THEN
58	Filtre_u(i,j)=0
59	ELSE
60	Filtre_u(i,j)=coeffu(index_vp(i),j)
61	ENDIF
62	ENDDO
63	ENDDO
64
65	DO j=1,jfiltnv
66	DO i=1,iim
67	IF (index_vp(i) < modfrstv(j)) THEN
68	Filtre_v(i,j)=0
69	ELSE
70	Filtre_v(i,j)=coeffv(index_vp(i),j)
71	ENDIF
72	ENDDO
73	ENDDO
74
75	DO j=jfiltsv,jjm
76	DO i=1,iim
77	IF (index_vp(i) < modfrstv(j)) THEN
78	Filtre_v(i,j)=0
79	ELSE
80	Filtre_v(i,j)=coeffv(index_vp(i),j)
81	ENDIF
82	ENDDO
83	ENDDO
84
85	DO j=2,jfiltnu
86	DO i=1,iim
87	IF (index_vp(i) < modfrstu(j)) THEN
88	Filtre_inv(i,j)=0
89	ELSE
90	Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
91	ENDIF
92	ENDDO
93	ENDDO
94
95	DO j=jfiltsu,jjm
96	DO i=1,iim
97	IF (index_vp(i) < modfrstu(j)) THEN
98	Filtre_inv(i,j)=0
99	ELSE
100	Filtre_inv(i,j)=coeffu(index_vp(i),j)/(1.+coeffu(index_vp(i),j))
101	ENDIF
102	ENDDO
103	ENDDO
104
105	#ifdef FFT_FFTW
106
107	WRITE (,)"COTH jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv"
108	WRITE (,)jfiltnu,jfiltsu,jfiltnv,jjm-jfiltsv
109	WRITE (,)MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1
110	CALL Init_FFT(iim,(llm+1)*(MAX(jfiltnu-2,jjm-jfiltsu,jfiltnv-2,jjm-jfiltsv)+1))
111	#else
112	CALL Init_FFT(iim,(jjm+1)*(llm+1))
113	#endif
114
115	END SUBROUTINE Init_filtre_fft
116
117	SUBROUTINE Filtre_u_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
118	USE mod_fft
119	#ifdef CPP_PARA
120	USE parallel_lmdz,ONLY : OMP_CHUNK
121	#endif
122	IMPLICIT NONE
123	include 'dimensions.h'
124	INTEGER,INTENT(IN) :: nlat
125	INTEGER,INTENT(IN) :: jj_begin
126	INTEGER,INTENT(IN) :: jj_end
127	INTEGER,INTENT(IN) :: nbniv
128	REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
129
130	REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
131	COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
132	INTEGER :: nb_vect
133	INTEGER :: i,j,l
134	INTEGER :: ll_nb
135
136	ll_nb=0
137	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
138	DO l=1,nbniv
139	ll_nb=ll_nb+1
140	DO j=1,jj_end-jj_begin+1
141	DO i=1,iim+1
142	vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
143	ENDDO
144	ENDDO
145	ENDDO
146	!$OMP END DO NOWAIT
147
148	nb_vect=(jj_end-jj_begin+1)*ll_nb
149
150	CALL FFT_forward(vect,TF_vect,nb_vect)
151
152	DO l=1,ll_nb
153	DO j=1,jj_end-jj_begin+1
154	DO i=1,iim/2+1
155	TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_u(i,jj_begin+j-1)
156	ENDDO
157	ENDDO
158	ENDDO
159
160	CALL FFT_backward(TF_vect,vect,nb_vect)
161
162
163	ll_nb=0
164	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
165	DO l=1,nbniv
166	ll_nb=ll_nb+1
167	DO j=1,jj_end-jj_begin+1
168	DO i=1,iim+1
169	vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
170	ENDDO
171	ENDDO
172	ENDDO
173	!$OMP END DO NOWAIT
174
175	END SUBROUTINE Filtre_u_fft
176
177
178	SUBROUTINE Filtre_v_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
179	USE mod_fft
180	#ifdef CPP_PARA
181	USE parallel_lmdz,ONLY : OMP_CHUNK
182	#endif
183	IMPLICIT NONE
184	INCLUDE 'dimensions.h'
185	INTEGER,INTENT(IN) :: nlat
186	INTEGER,INTENT(IN) :: jj_begin
187	INTEGER,INTENT(IN) :: jj_end
188	INTEGER,INTENT(IN) :: nbniv
189	REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
190
191	REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
192	COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
193	INTEGER :: nb_vect
194	INTEGER :: i,j,l
195	INTEGER :: ll_nb
196
197	ll_nb=0
198	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
199	DO l=1,nbniv
200	ll_nb=ll_nb+1
201	DO j=1,jj_end-jj_begin+1
202	DO i=1,iim+1
203	vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
204	ENDDO
205	ENDDO
206	ENDDO
207	!$OMP END DO NOWAIT
208
209
210	nb_vect=(jj_end-jj_begin+1)*ll_nb
211
212	CALL FFT_forward(vect,TF_vect,nb_vect)
213
214	DO l=1,ll_nb
215	DO j=1,jj_end-jj_begin+1
216	DO i=1,iim/2+1
217	TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_v(i,jj_begin+j-1)
218	ENDDO
219	ENDDO
220	ENDDO
221
222	CALL FFT_backward(TF_vect,vect,nb_vect)
223
224
225	ll_nb=0
226	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
227	DO l=1,nbniv
228	ll_nb=ll_nb+1
229	DO j=1,jj_end-jj_begin+1
230	DO i=1,iim+1
231	vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
232	ENDDO
233	ENDDO
234	ENDDO
235	!$OMP END DO NOWAIT
236
237	END SUBROUTINE Filtre_v_fft
238
239
240	SUBROUTINE Filtre_inv_fft(vect_inout,nlat,jj_begin,jj_end,nbniv)
241	USE mod_fft
242	#ifdef CPP_PARA
243	USE parallel_lmdz,ONLY : OMP_CHUNK
244	#endif
245	IMPLICIT NONE
246	INCLUDE 'dimensions.h'
247	INTEGER,INTENT(IN) :: nlat
248	INTEGER,INTENT(IN) :: jj_begin
249	INTEGER,INTENT(IN) :: jj_end
250	INTEGER,INTENT(IN) :: nbniv
251	REAL,INTENT(INOUT) :: vect_inout(iim+1,nlat,nbniv)
252
253	REAL :: vect(iim+inc,jj_end-jj_begin+1,nbniv)
254	COMPLEX :: TF_vect(iim/2+1,jj_end-jj_begin+1,nbniv)
255	INTEGER :: nb_vect
256	INTEGER :: i,j,l
257	INTEGER :: ll_nb
258
259	ll_nb=0
260	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
261	DO l=1,nbniv
262	ll_nb=ll_nb+1
263	DO j=1,jj_end-jj_begin+1
264	DO i=1,iim+1
265	vect(i,j,ll_nb)=vect_inout(i,j+jj_begin-1,l)
266	ENDDO
267	ENDDO
268	ENDDO
269	!$OMP END DO NOWAIT
270
271	nb_vect=(jj_end-jj_begin+1)*ll_nb
272
273	CALL FFT_forward(vect,TF_vect,nb_vect)
274
275	DO l=1,ll_nb
276	DO j=1,jj_end-jj_begin+1
277	DO i=1,iim/2+1
278	TF_vect(i,j,l)=TF_vect(i,j,l)*Filtre_inv(i,jj_begin+j-1)
279	ENDDO
280	ENDDO
281	ENDDO
282
283	CALL FFT_backward(TF_vect,vect,nb_vect)
284
285	ll_nb=0
286	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
287	DO l=1,nbniv
288	ll_nb=ll_nb+1
289	DO j=1,jj_end-jj_begin+1
290	DO i=1,iim+1
291	vect_inout(i,j+jj_begin-1,l)=vect(i,j,ll_nb)
292	ENDDO
293	ENDDO
294	ENDDO
295	!$OMP END DO NOWAIT
296
297	END SUBROUTINE Filtre_inv_fft
298
299	END MODULE mod_filtre_fft
300

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

Download in other formats:

Original Format