Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

filtreg.F @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago
Pluto PCM: Add the old Pluto LMDZ for reference (required prior step to making an LMDZ.PLUTO using the same framework as the other physics packages). TB+EM
Property svn:executable set to ``*
File size: 7.9 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	c write(,),'i,j,l',i,j,l
147	champ(i,j,l) = champ(i,j,l) * sdd1(i)
148	5 CONTINUE
149	c
150
151	IF( hemisph. EQ. 1 ) THEN
152
153	IF( ifiltre. EQ. -2 ) THEN
154	#ifdef CRAY
155	CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
156	* 1, iim, iim )
157	#else
158	#ifdef BLAS
159	CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
160	. champ(1,j,l), 1, 0.0, eignq, 1)
161	#else
162	DO k = 1, iim
163	eignq(k) = 0.0
164	ENDDO
165	DO k = 1, iim
166	DO i = 1, iim
167	eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l)
168	ENDDO
169	ENDDO
170	#endif
171	#endif
172	ELSE IF ( griscal ) THEN
173	#ifdef CRAY
174	CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
175	* 1, iim, iim )
176	#else
177	#ifdef BLAS
178	CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
179	. champ(1,j,l), 1, 0.0, eignq, 1)
180	#else
181	DO k = 1, iim
182	eignq(k) = 0.0
183	ENDDO
184	DO i = 1, iim
185	DO k = 1, iim
186	eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l)
187	ENDDO
188	ENDDO
189	#endif
190	#endif
191	ELSE
192	#ifdef CRAY
193	CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
194	* 1, iim, iim )
195	#else
196	#ifdef BLAS
197	CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
198	. champ(1,j,l), 1, 0.0, eignq, 1)
199	#else
200	DO k = 1, iim
201	eignq(k) = 0.0
202	ENDDO
203	DO i = 1, iim
204	DO k = 1, iim
205	eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l)
206	ENDDO
207	ENDDO
208	#endif
209	#endif
210	ENDIF
211
212	ELSE
213
214	IF( ifiltre. EQ. -2 ) THEN
215	#ifdef CRAY
216	CALL MXVA( matrinvs(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
217	* eignq, 1, iim, iim )
218	#else
219	#ifdef BLAS
220	CALL SGEMV("N", iim,iim, 1.0, matrinvs(1,1,j-jfiltsu+1),iim,
221	. champ(1,j,l), 1, 0.0, eignq, 1)
222	#else
223	DO k = 1, iim
224	eignq(k) = 0.0
225	ENDDO
226	DO i = 1, iim
227	DO k = 1, iim
228	eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l)
229	ENDDO
230	ENDDO
231	#endif
232	#endif
233	ELSE IF ( griscal ) THEN
234	#ifdef CRAY
235	CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
236	* eignq, 1, iim, iim )
237	#else
238	#ifdef BLAS
239	CALL SGEMV("N", iim,iim, 1.0, matriceus(1,1,j-jfiltsu+1),iim,
240	. champ(1,j,l), 1, 0.0, eignq, 1)
241	#else
242	DO k = 1, iim
243	eignq(k) = 0.0
244	ENDDO
245	DO i = 1, iim
246	DO k = 1, iim
247	eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l)
248	ENDDO
249	ENDDO
250	#endif
251	#endif
252	ELSE
253	#ifdef CRAY
254	CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 ,
255	* eignq, 1, iim, iim )
256	#else
257	#ifdef BLAS
258	CALL SGEMV("N", iim,iim, 1.0, matricevs(1,1,j-jfiltsv+1),iim,
259	. champ(1,j,l), 1, 0.0, eignq, 1)
260	#else
261	DO k = 1, iim
262	eignq(k) = 0.0
263	ENDDO
264	DO i = 1, iim
265	DO k = 1, iim
266	eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l)
267	ENDDO
268	ENDDO
269	#endif
270	#endif
271	ENDIF
272
273	ENDIF
274	c
275	IF( ifiltre.EQ. 2 ) THEN
276	DO 15 i = 1, iim
277	champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
278	15 CONTINUE
279	ELSE
280	DO 16 i=1,iim
281	champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
282	16 CONTINUE
283	ENDIF
284	c
285	champ( iip1,j,l ) = champ( 1,j,l )
286	c
287	30 CONTINUE
288	c
289	50 CONTINUE
290	c
291	100 CONTINUE
292	c
293	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
294	*filtrer, sur la grille des scalaires'/)
295	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
296	*ltrer, sur la grille de V ou de Z'/)
297	RETURN
298	END

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

Download in other formats:

Original Format