Context Navigation

filtreg.F @ 1086

Last change on this file since 1086 was 1086, checked in by yann meurdesoif, 15 years ago
Modifications Othman Bouzi : optimisation du filtre (remplacement opération matrice/vecteurs par matrice/matrice/matrice - BLAS), l'allocation des tableaux du filtre se fait maintenant dynamiquement (plus d'intervention manuelle dans parafiltre.h)
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 8.3 KB

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

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

Download in other formats:

Original Format