Context Navigation

filtreg.F @ 2

Last change on this file since 2 was 2, checked in by lmdz, 25 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 7.8 KB

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

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

Download in other formats:

Original Format