Context Navigation

filtreg.F @ 3403

Last change on this file since 3403 was 495, checked in by lmdzadmin, 21 years ago
Optimisation de differentes routines, IM, MAF, FH LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 10.7 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(2),jffil(2)
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	cIM REAL eignq(iim), sdd1(iim),sdd2(iim)
59	REAL eignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim)
60	LOGICAL griscal
61	INTEGER hemisph, iaire
62	c
63
64	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
65	* STOP'Pas de transformee simple dans cette version'
66
67	IF( iter.EQ. 2 ) THEN
68	PRINT *,' Pas d iteration du filtre dans cette version !'
69	* , ' Utiliser old_filtreg et repasser !'
70	STOP
71	ENDIF
72
73	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
74	PRINT *,' Cette routine ne calcule le filtre inverse que ',
75	* ' sur la grille des scalaires !'
76	STOP
77	ENDIF
78
79	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
80	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
81	*,' corriger et repasser !'
82	STOP
83	ENDIF
84	c
85
86	iim2 = iim * iim
87	immjm = iim * jjm
88	c
89	c
90	IF( griscal ) THEN
91	IF( nlat. NE. jjp1 ) THEN
92	PRINT 1111
93	STOP
94	ELSE
95	c
96	IF( iaire.EQ.1 ) THEN
97	CALL SCOPY( iim, sddv, 1, sdd1, 1 )
98	CALL SCOPY( iim, unsddv, 1, sdd2, 1 )
99	ELSE
100	CALL SCOPY( iim, unsddv, 1, sdd1, 1 )
101	CALL SCOPY( iim, sddv, 1, sdd2, 1 )
102	END IF
103	c
104	jdfil1 = 2
105	jffil1 = jfiltnu
106	jdfil2 = jfiltsu
107	jffil2 = jjm
108	END IF
109	ELSE
110	IF( nlat.NE.jjm ) THEN
111	PRINT 2222
112	STOP
113	ELSE
114	c
115	IF( iaire.EQ.1 ) THEN
116	CALL SCOPY( iim, sddu, 1, sdd1, 1 )
117	CALL SCOPY( iim, unsddu, 1, sdd2, 1 )
118	ELSE
119	CALL SCOPY( iim, unsddu, 1, sdd1, 1 )
120	CALL SCOPY( iim, sddu, 1, sdd2, 1 )
121	END IF
122	c
123	jdfil1 = 1
124	jffil1 = jfiltnv
125	jdfil2 = jfiltsv
126	jffil2 = jjm
127	END IF
128	END IF
129	c
130	jdfil(1) = jdfil1
131	jffil(1) = jffil1
132	jdfil(2) = jdfil2
133	jffil(2) = jffil2
134	c
135	DO 100 hemisph = 1, 2
136	c
137	c IF ( hemisph.EQ.1 ) THEN
138	c jdfil = jdfil1
139	c jffil = jffil1
140	c ELSE
141	c jdfil = jdfil2
142	c jffil = jffil2
143	c END IF
144
145
146	DO 50 l = 1, nbniv
147	DO 30 j = jdfil(hemisph),jffil(hemisph)
148
149
150	DO 5 i = 1, iim
151	champ(i,j,l) = champ(i,j,l) * sdd1(i)
152	5 CONTINUE
153	c
154	30 CONTINUE
155	50 CONTINUE
156
157	IF( hemisph. EQ. 1 ) THEN
158
159	IF( ifiltre. EQ. -2 ) THEN
160	#ifdef CRAY
161	DO l = 1, nbniv
162	DO j = jdfil(hemisph),jffil(hemisph)
163	CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
164	* 1, iim, iim )
165	ENDDO
166	ENDDO
167	#else
168	#ifdef BLAS
169	DO l = 1, nbniv
170	DO j = jdfil(hemisph),jffil(hemisph)
171	CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
172	. champ(1,j,l), 1, 0.0, eignq, 1)
173	ENDDO
174	ENDDO
175	#else
176	DO l = 1, nbniv
177	DO j = jdfil(hemisph),jffil(hemisph)
178	DO k = 1, iim
179	c eignq(k) = 0.0
180	eignq(k,j,l) = 0.0
181	ENDDO
182	DO k = 1, iim
183	DO i = 1, iim
184	c eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l)
185	eignq(k,j,l) = eignq(k,j,l) + matrinvn(k,i,j)*champ(i,j,l)
186	ENDDO
187	ENDDO
188	ENDDO
189	ENDDO
190	#endif
191	#endif
192	ELSE IF ( griscal ) THEN
193	#ifdef CRAY
194	DO l = 1, nbniv
195	DO j = jdfil(hemisph),jffil(hemisph)
196	CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
197	* 1, iim, iim )
198	ENDDO
199	ENDDO
200	#else
201	#ifdef BLAS
202	DO l = 1, nbniv
203	DO j = jdfil(hemisph),jffil(hemisph)
204	CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
205	. champ(1,j,l), 1, 0.0, eignq, 1)
206	ENDDO
207	ENDDO
208	#else
209	DO l = 1, nbniv
210	DO j = jdfil(hemisph),jffil(hemisph)
211	DO k = 1, iim
212	c eignq(k) = 0.0
213	eignq(k,j,l) = 0.0
214	ENDDO
215	DO i = 1, iim
216	DO k = 1, iim
217	c eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l)
218	eignq(k,j,l) = eignq(k,j,l) + matriceun(k,i,j)*champ(i,j,l)
219	ENDDO
220	ENDDO
221	ENDDO
222	ENDDO
223	#endif
224	#endif
225	ELSE
226	#ifdef CRAY
227	DO l = 1, nbniv
228	DO j = jdfil(hemisph),jffil(hemisph)
229	CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
230	* 1, iim, iim )
231	ENDDO
232	ENDDO
233	#else
234	#ifdef BLAS
235	DO l = 1, nbniv
236	DO j = jdfil(hemisph),jffil(hemisph)
237	CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
238	. champ(1,j,l), 1, 0.0, eignq, 1)
239	ENDDO
240	ENDDO
241	#else
242	DO l = 1, nbniv
243	DO j = jdfil(hemisph),jffil(hemisph)
244	DO k = 1, iim
245	c eignq(k) = 0.0
246	eignq(k,j,l) = 0.0
247	ENDDO
248	DO i = 1, iim
249	DO k = 1, iim
250	c eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l)
251	eignq(k,j,l) = eignq(k,j,l) + matricevn(k,i,j)*champ(i,j,l)
252	ENDDO
253	ENDDO
254	ENDDO
255	ENDDO
256	#endif
257	#endif
258	ENDIF
259
260	ELSE
261
262	IF( ifiltre. EQ. -2 ) THEN
263	#ifdef CRAY
264	DO l = 1, nbniv
265	DO j = jdfil(hemisph),jffil(hemisph)
266	CALL MXVA( matrinvs(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
267	* eignq, 1, iim, iim )
268	ENDDO
269	ENDDO
270	#else
271	#ifdef BLAS
272	DO l = 1, nbniv
273	DO j = jdfil(hemisph),jffil(hemisph)
274	CALL SGEMV("N", iim,iim, 1.0, matrinvs(1,1,j-jfiltsu+1),iim,
275	. champ(1,j,l), 1, 0.0, eignq, 1)
276	ENDDO
277	ENDDO
278	#else
279	DO l = 1, nbniv
280	DO j = jdfil(hemisph),jffil(hemisph)
281	DO k = 1, iim
282	c eignq(k) = 0.0
283	eignq(k,j,l) = 0.0
284	ENDDO
285	DO i = 1, iim
286	DO k = 1, iim
287	c eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l)
288	eignq(k,j,l) = eignq(k,j,l) +
289	.matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l)
290	ENDDO
291	ENDDO
292	ENDDO
293	ENDDO
294	#endif
295	#endif
296	ELSE IF ( griscal ) THEN
297	#ifdef CRAY
298	DO l = 1, nbniv
299	DO j = jdfil(hemisph),jffil(hemisph)
300	CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
301	* eignq, 1, iim, iim )
302	ENDDO
303	ENDDO
304	#else
305	#ifdef BLAS
306	DO l = 1, nbniv
307	DO j = jdfil(hemisph),jffil(hemisph)
308	CALL SGEMV("N", iim,iim, 1.0, matriceus(1,1,j-jfiltsu+1),iim,
309	. champ(1,j,l), 1, 0.0, eignq, 1)
310	ENDDO
311	ENDDO
312	#else
313	DO l = 1, nbniv
314	DO j = jdfil(hemisph),jffil(hemisph)
315	DO k = 1, iim
316	c eignq(k) = 0.0
317	eignq(k,j,l) = 0.0
318	ENDDO
319	DO i = 1, iim
320	DO k = 1, iim
321	c eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l)
322	eignq(k,j,l) = eignq(k,j,l) +
323	.matriceus(k,i,j-jfiltsu+1)*champ(i,j,l)
324	ENDDO
325	ENDDO
326	ENDDO
327	ENDDO
328	#endif
329	#endif
330	ELSE
331	#ifdef CRAY
332	DO l = 1, nbniv
333	DO j = jdfil(hemisph),jffil(hemisph)
334	CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 ,
335	* eignq, 1, iim, iim )
336	ENDDO
337	ENDDO
338	#else
339	#ifdef BLAS
340	DO l = 1, nbniv
341	DO j = jdfil(hemisph),jffil(hemisph)
342	CALL SGEMV("N", iim,iim, 1.0, matricevs(1,1,j-jfiltsv+1),iim,
343	. champ(1,j,l), 1, 0.0, eignq, 1)
344	ENDDO
345	ENDDO
346	#else
347	DO l = 1, nbniv
348	DO j = jdfil(hemisph),jffil(hemisph)
349	DO k = 1, iim
350	c eignq(k) = 0.0
351	eignq(k,j,l) = 0.0
352	ENDDO
353	DO i = 1, iim
354	DO k = 1, iim
355	c eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l)
356	eignq(k,j,l) = eignq(k,j,l) +
357	.matricevs(k,i,j-jfiltsv+1)*champ(i,j,l)
358	ENDDO
359	ENDDO
360	ENDDO
361	ENDDO
362	#endif
363	#endif
364	ENDIF
365
366	ENDIF
367	c
368	IF( ifiltre.EQ. 2 ) THEN
369	DO l = 1, nbniv
370	DO j = jdfil(hemisph),jffil(hemisph)
371	DO 15 i = 1, iim
372	c champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
373	champ( i,j,l ) = ( champ(i,j,l) + eignq(i,j,l) ) * sdd2(i)
374	15 CONTINUE
375	ENDDO
376	ENDDO
377	ELSE
378	DO l = 1, nbniv
379	DO j = jdfil(hemisph),jffil(hemisph)
380	DO 16 i=1,iim
381	c champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
382	champ( i,j,l ) = ( champ(i,j,l) - eignq(i,j,l) ) * sdd2(i)
383	16 CONTINUE
384	ENDDO
385	ENDDO
386	ENDIF
387	c
388	DO l = 1, nbniv
389	DO j = jdfil(hemisph),jffil(hemisph)
390	champ( iip1,j,l ) = champ( 1,j,l )
391	ENDDO
392	ENDDO
393	c
394	c 30 CONTINUE
395	c
396	c 50 CONTINUE
397	c
398	100 CONTINUE
399	c
400	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
401	*filtrer, sur la grille des scalaires'/)
402	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
403	*ltrer, sur la grille de V ou de Z'/)
404	RETURN
405	END

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

Download in other formats:

Original Format