Context Navigation

filtreg.F @ 2520

Last change on this file since 2520 was 704, checked in by Laurent Fairhead, 18 years ago
Inclusion des modifs de Y. Meurdesoif pour la version V3 LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 7.9 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE filtreg ( champ, nlat, nbniv, ifiltre,iaire,
5	. griscal ,iter)
6
7	IMPLICIT NONE
8	c=======================================================================
9	c
10	c Auteur: P. Le Van 07/10/97
11	c ------
12	c
13	c Objet: filtre matriciel longitudinal ,avec les matrices precalculees
14	c pour l'operateur Filtre .
15	c ------
16	c
17	c Arguments:
18	c ----------
19	c
20	c nblat nombre de latitudes a filtrer
21	c nbniv nombre de niveaux verticaux a filtrer
22	c champ(iip1,nblat,nbniv) en entree : champ a filtrer
23	c en sortie : champ filtre
24	c ifiltre +1 Transformee directe
25	c -1 Transformee inverse
26	c +2 Filtre directe
27	c -2 Filtre inverse
28	c
29	c iaire 1 si champ intensif
30	c 2 si champ extensif (pondere par les aires)
31	c
32	c iter 1 filtre simple
33	c
34	c=======================================================================
35	c
36	c
37	c Variable Intensive
38	c ifiltre = 1 filtre directe
39	c ifiltre =-1 filtre inverse
40	c
41	c Variable Extensive
42	c ifiltre = 2 filtre directe
43	c ifiltre =-2 filtre inverse
44	c
45	c
46	#include "dimensions.h"
47	#include "paramet.h"
48	#include "parafilt.h"
49	#include "coefils.h"
50	c
51	INTEGER nlat,nbniv,ifiltre,iter
52	INTEGER i,j,l,k
53	INTEGER iim2,immjm
54	INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
55
56	REAL champ( iip1,nlat,nbniv)
57	REAL matriceun,matriceus,matricevn,matricevs,matrinvn,matrinvs
58	COMMON/matrfil/matriceun(iim,iim,nfilun),matriceus(iim,iim,nfilus)
59	, , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs)
60	, , matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus)
61	REAL eignq(iim), sdd1(iim),sdd2(iim)
62	LOGICAL griscal
63	INTEGER hemisph, iaire
64	c
65
66	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
67	* STOP'Pas de transformee simple dans cette version'
68
69	IF( iter.EQ. 2 ) THEN
70	PRINT *,' Pas d iteration du filtre dans cette version !'
71	* , ' Utiliser old_filtreg et repasser !'
72	STOP
73	ENDIF
74
75	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
76	PRINT *,' Cette routine ne calcule le filtre inverse que ',
77	* ' sur la grille des scalaires !'
78	STOP
79	ENDIF
80
81	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
82	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
83	*,' corriger et repasser !'
84	STOP
85	ENDIF
86	c
87
88	iim2 = iim * iim
89	immjm = iim * jjm
90	c
91	c
92	IF( griscal ) THEN
93	IF( nlat. NE. jjp1 ) THEN
94	PRINT 1111
95	STOP
96	ELSE
97	c
98	IF( iaire.EQ.1 ) THEN
99	CALL SCOPY( iim, sddv, 1, sdd1, 1 )
100	CALL SCOPY( iim, unsddv, 1, sdd2, 1 )
101	ELSE
102	CALL SCOPY( iim, unsddv, 1, sdd1, 1 )
103	CALL SCOPY( iim, sddv, 1, sdd2, 1 )
104	END IF
105	c
106	jdfil1 = 2
107	jffil1 = jfiltnu
108	jdfil2 = jfiltsu
109	jffil2 = jjm
110	END IF
111	ELSE
112	IF( nlat.NE.jjm ) THEN
113	PRINT 2222
114	STOP
115	ELSE
116	c
117	IF( iaire.EQ.1 ) THEN
118	CALL SCOPY( iim, sddu, 1, sdd1, 1 )
119	CALL SCOPY( iim, unsddu, 1, sdd2, 1 )
120	ELSE
121	CALL SCOPY( iim, unsddu, 1, sdd1, 1 )
122	CALL SCOPY( iim, sddu, 1, sdd2, 1 )
123	END IF
124	c
125	jdfil1 = 1
126	jffil1 = jfiltnv
127	jdfil2 = jfiltsv
128	jffil2 = jjm
129	END IF
130	END IF
131	c
132	c
133	DO 100 hemisph = 1, 2
134	c
135	IF ( hemisph.EQ.1 ) THEN
136	jdfil = jdfil1
137	jffil = jffil1
138	ELSE
139	jdfil = jdfil2
140	jffil = jffil2
141	END IF
142
143
144	DO 50 l = 1, nbniv
145	DO 30 j = jdfil,jffil
146
147
148	DO 5 i = 1, iim
149	champ(i,j,l) = champ(i,j,l) * sdd1(i)
150	5 CONTINUE
151	c
152
153	IF( hemisph. EQ. 1 ) THEN
154
155	IF( ifiltre. EQ. -2 ) THEN
156	#ifdef CRAY
157	CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
158	* 1, iim, iim )
159	#else
160	#ifdef BLAS
161	CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
162	. champ(1,j,l), 1, 0.0, eignq, 1)
163	#else
164	DO k = 1, iim
165	eignq(k) = 0.0
166	ENDDO
167	DO k = 1, iim
168	DO i = 1, iim
169	eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l)
170	ENDDO
171	ENDDO
172	#endif
173	#endif
174	ELSE IF ( griscal ) THEN
175	#ifdef CRAY
176	CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
177	* 1, iim, iim )
178	#else
179	#ifdef BLAS
180	CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
181	. champ(1,j,l), 1, 0.0, eignq, 1)
182	#else
183	DO k = 1, iim
184	eignq(k) = 0.0
185	ENDDO
186	DO i = 1, iim
187	DO k = 1, iim
188	eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l)
189	ENDDO
190	ENDDO
191	#endif
192	#endif
193	ELSE
194	#ifdef CRAY
195	CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
196	* 1, iim, iim )
197	#else
198	#ifdef BLAS
199	CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
200	. champ(1,j,l), 1, 0.0, eignq, 1)
201	#else
202	DO k = 1, iim
203	eignq(k) = 0.0
204	ENDDO
205	DO i = 1, iim
206	DO k = 1, iim
207	eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l)
208	ENDDO
209	ENDDO
210	#endif
211	#endif
212	ENDIF
213
214	ELSE
215
216	IF( ifiltre. EQ. -2 ) THEN
217	#ifdef CRAY
218	CALL MXVA( matrinvs(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
219	* eignq, 1, iim, iim )
220	#else
221	#ifdef BLAS
222	CALL SGEMV("N", iim,iim, 1.0, matrinvs(1,1,j-jfiltsu+1),iim,
223	. champ(1,j,l), 1, 0.0, eignq, 1)
224	#else
225	DO k = 1, iim
226	eignq(k) = 0.0
227	ENDDO
228	DO i = 1, iim
229	DO k = 1, iim
230	eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l)
231	ENDDO
232	ENDDO
233	#endif
234	#endif
235	ELSE IF ( griscal ) THEN
236	#ifdef CRAY
237	CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
238	* eignq, 1, iim, iim )
239	#else
240	#ifdef BLAS
241	CALL SGEMV("N", iim,iim, 1.0, matriceus(1,1,j-jfiltsu+1),iim,
242	. champ(1,j,l), 1, 0.0, eignq, 1)
243	#else
244	DO k = 1, iim
245	eignq(k) = 0.0
246	ENDDO
247	DO i = 1, iim
248	DO k = 1, iim
249	eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l)
250	ENDDO
251	ENDDO
252	#endif
253	#endif
254	ELSE
255	#ifdef CRAY
256	CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 ,
257	* eignq, 1, iim, iim )
258	#else
259	#ifdef BLAS
260	CALL SGEMV("N", iim,iim, 1.0, matricevs(1,1,j-jfiltsv+1),iim,
261	. champ(1,j,l), 1, 0.0, eignq, 1)
262	#else
263	DO k = 1, iim
264	eignq(k) = 0.0
265	ENDDO
266	DO i = 1, iim
267	DO k = 1, iim
268	eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l)
269	ENDDO
270	ENDDO
271	#endif
272	#endif
273	ENDIF
274
275	ENDIF
276	c
277	IF( ifiltre.EQ. 2 ) THEN
278	DO 15 i = 1, iim
279	champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
280	15 CONTINUE
281	ELSE
282	DO 16 i=1,iim
283	champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
284	16 CONTINUE
285	ENDIF
286	c
287	champ( iip1,j,l ) = champ( 1,j,l )
288	c
289	30 CONTINUE
290	c
291	50 CONTINUE
292	c
293	100 CONTINUE
294	c
295	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
296	*filtrer, sur la grille des scalaires'/)
297	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
298	*ltrer, sur la grille de V ou de Z'/)
299	RETURN
300	END

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

Download in other formats:

Original Format