Context Navigation

filtreg_p.F @ 1068

Last change on this file since 1068 was 985, checked in by Laurent Fairhead, 16 years ago
Mise a jour de dyn3dpar par rapport a dyn3d, inclusion OpenMP et filtre FFT YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 10.0 KB

Line
1
2
3	SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv,
4	. ifiltre, iaire, griscal ,iter)
5	USE Parallel, only : OMP_CHUNK
6	USE mod_filtre_fft
7	USE timer_filtre
8	IMPLICIT NONE
9
10	c=======================================================================
11	c
12	c Auteur: P. Le Van 07/10/97
13	c ------
14	c
15	c Objet: filtre matriciel longitudinal ,avec les matrices precalculees
16	c pour l'operateur Filtre .
17	c ------
18	c
19	c Arguments:
20	c ----------
21	c
22	c
23	c ibeg..iend lattitude a filtrer
24	c nlat nombre de latitudes du champ
25	c nbniv nombre de niveaux verticaux a filtrer
26	c champ(iip1,nblat,nbniv) en entree : champ a filtrer
27	c en sortie : champ filtre
28	c ifiltre +1 Transformee directe
29	c -1 Transformee inverse
30	c +2 Filtre directe
31	c -2 Filtre inverse
32	c
33	c iaire 1 si champ intensif
34	c 2 si champ extensif (pondere par les aires)
35	c
36	c iter 1 filtre simple
37	c
38	c=======================================================================
39	c
40	c
41	c Variable Intensive
42	c ifiltre = 1 filtre directe
43	c ifiltre =-1 filtre inverse
44	c
45	c Variable Extensive
46	c ifiltre = 2 filtre directe
47	c ifiltre =-2 filtre inverse
48	c
49	c
50	#include "dimensions.h"
51	#include "paramet.h"
52	#include "parafilt.h"
53	#include "coefils.h"
54	c
55	INTEGER ibeg,iend,nlat,nbniv,ifiltre,iter
56	INTEGER i,j,l,k
57	INTEGER iim2,immjm
58	INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
59
60	REAL champ( iip1,nlat,nbniv)
61	REAL matriceun,matriceus,matricevn,matricevs,matrinvn,matrinvs
62	COMMON/matrfil/matriceun(iim,iim,nfilun),matriceus(iim,iim,nfilus)
63	, , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs)
64	, , matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus)
65	cym REAL eignq(iim), sdd1(iim),sdd2(iim)
66
67	REAL eignq(iim)
68	REAL :: sdd1(iim),sdd2(iim)
69
70	LOGICAL griscal
71	INTEGER hemisph, iaire
72
73	REAL :: champ_fft(iip1,nlat,nbniv)
74	REAL :: champ_in(iip1,nlat,nbniv)
75
76	REAL,SAVE,TARGET :: sddu_loc(iim)
77	REAL,SAVE,TARGET :: sddv_loc(iim)
78	REAL,SAVE,TARGET :: unsddu_loc(iim)
79	REAL,SAVE,TARGET :: unsddv_loc(iim)
80	c$OMP THREADPRIVATE(sddu_loc,sddv_loc,unsddu_loc,unsddv_loc)
81	LOGICAL,SAVE :: first=.TRUE.
82	c$OMP THREADPRIVATE(first)
83
84	IF (first) THEN
85	sddu_loc(1:iim)=sddu(1:iim)
86	sddv_loc(1:iim)=sddv(1:iim)
87	unsddu_loc(1:iim)=unsddu(1:iim)
88	unsddv_loc(1:iim)=unsddv(1:iim)
89	CALL Init_timer
90	first=.FALSE.
91	c PRINT *,"----> sddu_loc=",sddu_loc
92	c PRINT *,"----> sddv_loc=",sddv_loc
93	c PRINT *,"----> unsddu_loc=",unsddu_loc
94	c PRINT *,"----> unsddv_loc=",unsddv_loc
95	ENDIF
96
97	c$OMP MASTER
98	CALL start_timer
99	c$OMP END MASTER
100
101	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
102	* STOP'Pas de transformee simple dans cette version'
103
104	IF( iter.EQ. 2 ) THEN
105	PRINT *,' Pas d iteration du filtre dans cette version !'
106	* , ' Utiliser old_filtreg et repasser !'
107	STOP
108	ENDIF
109
110	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
111	PRINT *,' Cette routine ne calcule le filtre inverse que ',
112	* ' sur la grille des scalaires !'
113	STOP
114	ENDIF
115
116	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
117	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
118	*,' corriger et repasser !'
119	STOP
120	ENDIF
121	c
122
123	iim2 = iim * iim
124	immjm = iim * jjm
125	c
126	c
127	IF( griscal ) THEN
128	IF( nlat. NE. jjp1 ) THEN
129	PRINT 1111
130	STOP
131	ELSE
132	c
133	IF( iaire.EQ.1 ) THEN
134	cym CALL SCOPY( iim, sddv, 1, sdd1, 1 )
135	cym CALL SCOPY( iim, unsddv, 1, sdd2, 1 )
136	cym sdd1=>sddv_loc
137	cym sdd2=>unsddv_loc
138	sdd1(1:iim)=sddv_loc(1:iim)
139	sdd2(1:iim)=unsddv_loc(1:iim)
140	ELSE
141	cym CALL SCOPY( iim, unsddv, 1, sdd1, 1 )
142	cym CALL SCOPY( iim, sddv, 1, sdd2, 1 )
143	sdd1(1:iim)=unsddv_loc(1:iim)
144	sdd2(1:iim)=sddv_loc(1:iim)
145	END IF
146	c
147	jdfil1 = 2
148	jffil1 = jfiltnu
149	jdfil2 = jfiltsu
150	jffil2 = jjm
151	END IF
152	ELSE
153	IF( nlat.NE.jjm ) THEN
154	PRINT 2222
155	STOP
156	ELSE
157	c
158	IF( iaire.EQ.1 ) THEN
159	cym CALL SCOPY( iim, sddu, 1, sdd1, 1 )
160	cym CALL SCOPY( iim, unsddu, 1, sdd2, 1 )
161	cym sdd1=>sddu_loc
162	cym sdd2=>unsddu_loc
163	sdd1(1:iim)=sddu_loc(1:iim)
164	sdd2(1:iim)=unsddu_loc(1:iim)
165
166	ELSE
167	cym CALL SCOPY( iim, unsddu, 1, sdd1, 1 )
168	cym CALL SCOPY( iim, sddu, 1, sdd2, 1 )
169	cym sdd1=>unsddu_loc
170	cym sdd2=>sddu_loc
171	sdd1(1:iim)=unsddu_loc(1:iim)
172	sdd2(1:iim)=sddu_loc(1:iim)
173	END IF
174	c
175	jdfil1 = 1
176	jffil1 = jfiltnv
177	jdfil2 = jfiltsv
178	jffil2 = jjm
179	END IF
180	END IF
181
182	c PRINT *,"APPEL a filtreg --> sdd1=",sdd1
183	c PRINT *,"APPEL a filtreg --> sdd2=",sdd2
184	c PRINT *,"----> sddu_loc=",sddu_loc
185	c PRINT *,"----> sddv_loc=",sddv_loc
186	c PRINT *,"----> unsddu_loc=",unsddu_loc
187	c PRINT *,"----> unsddv_loc=",unsddv_loc
188
189	c
190	c
191	DO 100 hemisph = 1, 2
192	c
193	IF ( hemisph.EQ.1 ) THEN
194	c ym
195	jdfil = max(jdfil1,ibeg)
196	jffil = min(jffil1,iend)
197	ELSE
198	c ym
199	jdfil = max(jdfil2,ibeg)
200	jffil = min(jffil2,iend)
201	END IF
202
203
204	cccccccccccccccccccccccccccccccccccccccccccc
205	c Utilisation du filtre classique
206	cccccccccccccccccccccccccccccccccccccccccccc
207
208	IF (.NOT. use_filtre_fft) THEN
209
210	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
211	DO 50 l = 1, nbniv
212	DO 30 j = jdfil,jffil
213
214
215	DO 5 i = 1, iim
216	champ(i,j,l) = champ(i,j,l) * sdd1(i)
217	5 CONTINUE
218	c
219
220	IF( hemisph. EQ. 1 ) THEN
221
222	IF( ifiltre. EQ. -2 ) THEN
223
224
225	CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
226	. champ(1,j,l), 1, 0.0, eignq, 1)
227
228
229	ELSE IF ( griscal ) THEN
230
231	CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
232	. champ(1,j,l), 1, 0.0, eignq, 1)
233
234	ELSE
235
236	CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
237	. champ(1,j,l), 1, 0.0, eignq, 1)
238	ENDIF
239
240	ELSE
241
242	IF( ifiltre. EQ. -2 ) THEN
243
244	CALL SGEMV("N",iim,iim,1.0, matrinvs(1,1,j-jfiltsu+1),iim,
245	. champ(1,j,l), 1, 0.0, eignq, 1)
246
247	ELSE IF ( griscal ) THEN
248
249	CALL SGEMV("N",iim,iim,1.0,matriceus(1,1,j-jfiltsu+1),iim,
250	. champ(1,j,l), 1, 0.0, eignq, 1)
251	ELSE
252
253	CALL SGEMV("N",iim,iim,1.0,matricevs(1,1,j-jfiltsv+1),iim,
254	. champ(1,j,l), 1, 0.0, eignq, 1)
255	ENDIF
256
257	ENDIF
258
259
260	c
261	IF( ifiltre.EQ. 2 ) THEN
262
263	DO 15 i = 1, iim
264	champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
265	15 CONTINUE
266
267	ELSE
268
269	DO 16 i=1,iim
270	champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
271	16 CONTINUE
272
273	ENDIF
274	c
275	champ( iip1,j,l ) = champ( 1,j,l )
276	c
277	30 CONTINUE
278	c
279	50 CONTINUE
280	c$OMP END DO NOWAIT
281
282	ccccccccccccccccccccccccccccccccccccccccccccc
283	c Utilisation du filtre FFT
284	ccccccccccccccccccccccccccccccccccccccccccccc
285
286	ELSE
287
288	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
289	DO l=1,nbniv
290	DO j=jdfil,jffil
291	DO i = 1, iim
292	champ( i,j,l)= champ(i,j,l)*sdd1(i)
293	champ_fft( i,j,l) = champ(i,j,l)
294	ENDDO
295	ENDDO
296	ENDDO
297	c$OMP END DO NOWAIT
298
299	IF (jdfil<=jffil) THEN
300	IF( ifiltre. EQ. -2 ) THEN
301	CALL Filtre_inv_fft(champ_fft,nlat,jdfil,jffil,nbniv)
302	ELSE IF ( griscal ) THEN
303	CALL Filtre_u_fft(champ_fft,nlat,jdfil,jffil,nbniv)
304	ELSE
305	CALL Filtre_v_fft(champ_fft,nlat,jdfil,jffil,nbniv)
306	ENDIF
307	ENDIF
308
309
310	IF( ifiltre.EQ. 2 ) THEN
311	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
312	DO l=1,nbniv
313	DO j=jdfil,jffil
314	DO i = 1, iim
315	champ( i,j,l)=(champ(i,j,l)+champ_fft(i,j,l))
316	. *sdd2(i)
317	ENDDO
318	ENDDO
319	ENDDO
320	c$OMP END DO NOWAIT
321	ELSE
322
323	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
324	DO l=1,nbniv
325	DO j=jdfil,jffil
326	DO i = 1, iim
327	champ(i,j,l)=(champ(i,j,l)-champ_fft(i,j,l))
328	. *sdd2(i)
329	ENDDO
330	ENDDO
331	ENDDO
332	c$OMP END DO NOWAIT
333	ENDIF
334	c
335	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
336	DO l=1,nbniv
337	DO j=jdfil,jffil
338	! champ_FFT( iip1,j,l ) = champ_FFT( 1,j,l )
339	champ( iip1,j,l ) = champ( 1,j,l )
340	ENDDO
341	ENDDO
342	c$OMP END DO NOWAIT
343	ENDIF
344	c Fin de la zone de filtrage
345
346
347	100 CONTINUE
348
349	! DO j=1,nlat
350	!
351	! PRINT *,"check FFT ----> Delta(",j,")=",
352	! & sum(champ(:,j,:)-champ_fft(:,j,:))/sum(champ(:,j,:)),
353	! & sum(champ(:,j,:)-champ_in(:,j,:))/sum(champ(:,j,:))
354	! ENDDO
355
356	! PRINT *,"check FFT ----> Delta(",j,")=",
357	! & sum(champ-champ_fft)/sum(champ)
358	!
359
360	c
361	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
362	*filtrer, sur la grille des scalaires'/)
363	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
364	*ltrer, sur la grille de V ou de Z'/)
365	c$OMP MASTER
366	CALL stop_timer
367	c$OMP END MASTER
368	RETURN
369	END

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

Download in other formats:

Original Format