Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

filtreg_p.F @ 5415

Last change on this file since 5415 was 708, checked in by Laurent Fairhead, 18 years ago
Versions parallèlisées des routines YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 8.1 KB

Line
1	SUBROUTINE filtreg_p ( champ, ibeg, iend, nlat, nbniv,
2	. ifiltre, iaire, griscal ,iter)
3	USE Parallel, only : OMP_CHUNK
4	IMPLICIT NONE
5
6	c=======================================================================
7	c
8	c Auteur: P. Le Van 07/10/97
9	c ------
10	c
11	c Objet: filtre matriciel longitudinal ,avec les matrices precalculees
12	c pour l'operateur Filtre .
13	c ------
14	c
15	c Arguments:
16	c ----------
17	c
18	c
19	c ibeg..iend lattitude a filtrer
20	c nlat nombre de latitudes du champ
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 ibeg,iend,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	c ym
137	jdfil = max(jdfil1,ibeg)
138	jffil = min(jffil1,iend)
139	ELSE
140	c ym
141	jdfil = max(jdfil2,ibeg)
142	jffil = min(jffil2,iend)
143	END IF
144
145	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
146	DO 50 l = 1, nbniv
147	DO 30 j = jdfil,jffil
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
155	IF( hemisph. EQ. 1 ) THEN
156
157	IF( ifiltre. EQ. -2 ) THEN
158	#ifdef CRAY
159	CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
160	* 1, iim, iim )
161	#else
162	#ifdef BLAS
163	CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
164	. champ(1,j,l), 1, 0.0, eignq, 1)
165	#else
166	DO k = 1, iim
167	eignq(k) = 0.0
168	ENDDO
169	DO k = 1, iim
170	DO i = 1, iim
171	eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l)
172	ENDDO
173	ENDDO
174	#endif
175	#endif
176	ELSE IF ( griscal ) THEN
177	#ifdef CRAY
178	CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
179	* 1, iim, iim )
180	#else
181	#ifdef BLAS
182	CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
183	. champ(1,j,l), 1, 0.0, eignq, 1)
184	#else
185	DO k = 1, iim
186	eignq(k) = 0.0
187	ENDDO
188	DO i = 1, iim
189	DO k = 1, iim
190	eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l)
191	ENDDO
192	ENDDO
193	#endif
194	#endif
195	ELSE
196	#ifdef CRAY
197	CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
198	* 1, iim, iim )
199	#else
200	#ifdef BLAS
201	CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
202	. champ(1,j,l), 1, 0.0, eignq, 1)
203	#else
204	DO k = 1, iim
205	eignq(k) = 0.0
206	ENDDO
207	DO i = 1, iim
208	DO k = 1, iim
209	eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l)
210	ENDDO
211	ENDDO
212	#endif
213	#endif
214	ENDIF
215
216	ELSE
217
218	IF( ifiltre. EQ. -2 ) THEN
219	#ifdef CRAY
220	CALL MXVA( matrinvs(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
221	* eignq, 1, iim, iim )
222	#else
223	#ifdef BLAS
224	CALL SGEMV("N", iim,iim, 1.0, matrinvs(1,1,j-jfiltsu+1),iim,
225	. champ(1,j,l), 1, 0.0, eignq, 1)
226	#else
227	DO k = 1, iim
228	eignq(k) = 0.0
229	ENDDO
230	DO i = 1, iim
231	DO k = 1, iim
232	eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l)
233	ENDDO
234	ENDDO
235	#endif
236	#endif
237	ELSE IF ( griscal ) THEN
238	#ifdef CRAY
239	CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
240	* eignq, 1, iim, iim )
241	#else
242	#ifdef BLAS
243	CALL SGEMV("N", iim,iim, 1.0, matriceus(1,1,j-jfiltsu+1),iim,
244	. champ(1,j,l), 1, 0.0, eignq, 1)
245	#else
246	DO k = 1, iim
247	eignq(k) = 0.0
248	ENDDO
249	DO i = 1, iim
250	DO k = 1, iim
251	eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l)
252	ENDDO
253	ENDDO
254	#endif
255	#endif
256	ELSE
257	#ifdef CRAY
258	CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 ,
259	* eignq, 1, iim, iim )
260	#else
261	#ifdef BLAS
262	CALL SGEMV("N", iim,iim, 1.0, matricevs(1,1,j-jfiltsv+1),iim,
263	. champ(1,j,l), 1, 0.0, eignq, 1)
264	#else
265	DO k = 1, iim
266	eignq(k) = 0.0
267	ENDDO
268	DO i = 1, iim
269	DO k = 1, iim
270	eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l)
271	ENDDO
272	ENDDO
273	#endif
274	#endif
275	ENDIF
276
277	ENDIF
278	c
279	IF( ifiltre.EQ. 2 ) THEN
280	DO 15 i = 1, iim
281	champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
282	15 CONTINUE
283	ELSE
284	DO 16 i=1,iim
285	champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
286	16 CONTINUE
287	ENDIF
288	c
289	champ( iip1,j,l ) = champ( 1,j,l )
290	c
291	30 CONTINUE
292	c
293	50 CONTINUE
294	c$OMP END DO NOWAIT
295	c
296	100 CONTINUE
297	c
298	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
299	*filtrer, sur la grille des scalaires'/)
300	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
301	*ltrer, sur la grille de V ou de Z'/)
302	RETURN
303	END

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

Download in other formats:

Original Format