Context Navigation

filtreg.F @ 797

Last change on this file since 797 was 524, checked in by lmdzadmin, 21 years ago
Initial revision
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