Context Navigation

filtreg.F @ 5456

Last change on this file since 5456 was 4440, checked in by fhourdin, 23 months ago
Deactivating filtre and check if iim=1 (2D y-z)
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 9.2 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE filtreg ( champ, nlat, nbniv, ifiltre,iaire,
5	& griscal ,iter)
6
7	USE filtreg_mod
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 ( iim == 1 ) return ! no filtre in 2D y-z
75
76	IF (first) THEN
77	sdd12(1:iim,type_sddu) = sddu(1:iim)
78	sdd12(1:iim,type_sddv) = sddv(1:iim)
79	sdd12(1:iim,type_unsddu) = unsddu(1:iim)
80	sdd12(1:iim,type_unsddv) = unsddv(1:iim)
81
82	first=.FALSE.
83	ENDIF
84
85	IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
86	& STOP'Pas de transformee simple dans cette version'
87
88	IF( iter.EQ. 2 ) THEN
89	PRINT *,' Pas d iteration du filtre dans cette version !'
90	& , ' Utiliser old_filtreg et repasser !'
91	STOP
92	ENDIF
93
94	IF( ifiltre.EQ. -2 .AND..NOT.griscal ) THEN
95	PRINT *,' Cette routine ne calcule le filtre inverse que '
96	& , ' sur la grille des scalaires !'
97	STOP
98	ENDIF
99
100	IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 ) THEN
101	PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
102	& , ' corriger et repasser !'
103	STOP
104	ENDIF
105
106	iim2 = iim * iim
107	immjm = iim * jjm
108
109	IF( griscal ) THEN
110	IF( nlat. NE. jjp1 ) THEN
111	PRINT 1111
112	STOP
113	ELSE
114
115	IF( iaire.EQ.1 ) THEN
116	sdd1_type = type_sddv
117	sdd2_type = type_unsddv
118	ELSE
119	sdd1_type = type_unsddv
120	sdd2_type = type_sddv
121	ENDIF
122
123	c IF( iaire.EQ.1 ) THEN
124	c CALL SCOPY( iim, sddv, 1, sdd1, 1 )
125	c CALL SCOPY( iim, unsddv, 1, sdd2, 1 )
126	c ELSE
127	c CALL SCOPY( iim, unsddv, 1, sdd1, 1 )
128	c CALL SCOPY( iim, sddv, 1, sdd2, 1 )
129	c END IF
130
131	jdfil1 = 2
132	jffil1 = jfiltnu
133	jdfil2 = jfiltsu
134	jffil2 = jjm
135	END IF
136	ELSE
137	IF( nlat.NE.jjm ) THEN
138	PRINT 2222
139	STOP
140	ELSE
141
142	IF( iaire.EQ.1 ) THEN
143	sdd1_type = type_sddu
144	sdd2_type = type_unsddu
145	ELSE
146	sdd1_type = type_unsddu
147	sdd2_type = type_sddu
148	ENDIF
149
150	c IF( iaire.EQ.1 ) THEN
151	c CALL SCOPY( iim, sddu, 1, sdd1, 1 )
152	c CALL SCOPY( iim, unsddu, 1, sdd2, 1 )
153	c ELSE
154	c CALL SCOPY( iim, unsddu, 1, sdd1, 1 )
155	c CALL SCOPY( iim, sddu, 1, sdd2, 1 )
156	c END IF
157
158	jdfil1 = 1
159	jffil1 = jfiltnv
160	jdfil2 = jfiltsv
161	jffil2 = jjm
162	END IF
163	END IF
164
165	DO hemisph = 1, 2
166
167	IF ( hemisph.EQ.1 ) THEN
168	jdfil = jdfil1
169	jffil = jffil1
170	ELSE
171	jdfil = jdfil2
172	jffil = jffil2
173	END IF
174
175	DO l = 1, nbniv
176	DO j = jdfil,jffil
177	DO i = 1, iim
178	champ(i,j,l) = champ(i,j,l) * sdd12(i,sdd1_type) ! sdd1(i)
179	END DO
180	END DO
181	END DO
182
183	IF( hemisph. EQ. 1 ) THEN
184
185	IF( ifiltre. EQ. -2 ) THEN
186
187	DO j = jdfil,jffil
188	#ifdef BLAS
189	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
190	& matrinvn(1,1,j),
191	& iim, champ(1,j,1), iip1*nlat, 0.0,
192	& eignq(1,j-jdfil+1,1), iim*nlat)
193	#else
194	eignq(:,j-jdfil+1,:)
195	$ = matmul(matrinvn(:,:,j), champ(:iim,j,:))
196	#endif
197	END DO
198
199	ELSE IF ( griscal ) THEN
200
201	DO j = jdfil,jffil
202	#ifdef BLAS
203	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
204	& matriceun(1,1,j),
205	& iim, champ(1,j,1), iip1*nlat, 0.0,
206	& eignq(1,j-jdfil+1,1), iim*nlat)
207	#else
208	eignq(:,j-jdfil+1,:)
209	$ = matmul(matriceun(:,:,j), champ(:iim,j,:))
210	#endif
211	END DO
212
213	ELSE
214
215	DO j = jdfil,jffil
216	#ifdef BLAS
217	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
218	& matricevn(1,1,j),
219	& iim, champ(1,j,1), iip1*nlat, 0.0,
220	& eignq(1,j-jdfil+1,1), iim*nlat)
221	#else
222	eignq(:,j-jdfil+1,:)
223	$ = matmul(matricevn(:,:,j), champ(:iim,j,:))
224	#endif
225	END DO
226
227	ENDIF
228
229	ELSE
230
231	IF( ifiltre. EQ. -2 ) THEN
232
233	DO j = jdfil,jffil
234	#ifdef BLAS
235	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
236	& matrinvs(1,1,j-jfiltsu+1),
237	& iim, champ(1,j,1), iip1*nlat, 0.0,
238	& eignq(1,j-jdfil+1,1), iim*nlat)
239	#else
240	eignq(:,j-jdfil+1,:)
241	$ = matmul(matrinvs(:,:,j-jfiltsu+1),
242	$ champ(:iim,j,:))
243	#endif
244	END DO
245
246
247	ELSE IF ( griscal ) THEN
248
249	DO j = jdfil,jffil
250	#ifdef BLAS
251	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
252	& matriceus(1,1,j-jfiltsu+1),
253	& iim, champ(1,j,1), iip1*nlat, 0.0,
254	& eignq(1,j-jdfil+1,1), iim*nlat)
255	#else
256	eignq(:,j-jdfil+1,:)
257	$ = matmul(matriceus(:,:,j-jfiltsu+1),
258	$ champ(:iim,j,:))
259	#endif
260	END DO
261
262	ELSE
263
264	DO j = jdfil,jffil
265	#ifdef BLAS
266	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
267	& matricevs(1,1,j-jfiltsv+1),
268	& iim, champ(1,j,1), iip1*nlat, 0.0,
269	& eignq(1,j-jdfil+1,1), iim*nlat)
270	#else
271	eignq(:,j-jdfil+1,:)
272	$ = matmul(matricevs(:,:,j-jfiltsv+1),
273	$ champ(:iim,j,:))
274	#endif
275	END DO
276
277	ENDIF
278
279	ENDIF
280
281	IF( ifiltre.EQ. 2 ) THEN
282
283	DO l = 1, nbniv
284	DO j = jdfil,jffil
285	DO i = 1, iim
286	champ( i,j,l ) =
287	& (champ(i,j,l) + eignq(i,j-jdfil+1,l))
288	& * sdd12(i,sdd2_type) ! sdd2(i)
289	END DO
290	END DO
291	END DO
292
293	ELSE
294
295	DO l = 1, nbniv
296	DO j = jdfil,jffil
297	DO i = 1, iim
298	champ( i,j,l ) =
299	& (champ(i,j,l) - eignq(i,j-jdfil+1,l))
300	& * sdd12(i,sdd2_type) ! sdd2(i)
301	END DO
302	END DO
303	END DO
304
305	ENDIF
306
307	DO l = 1, nbniv
308	DO j = jdfil,jffil
309	champ( iip1,j,l ) = champ( 1,j,l )
310	END DO
311	END DO
312
313
314	ENDDO
315
316	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
317	& filtrer, sur la grille des scalaires'/)
318	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
319	& ltrer, sur la grille de V ou de Z'/)
320	RETURN
321	END

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

Download in other formats:

Original Format