Context Navigation

filtreg.F @ 2780

Last change on this file since 2780 was 2626, checked in by emillour, 4 years ago
Common dynamics: Some rather harmless OpenMP fixes in the filtering identified by recent versions of ifort. EM
File size: 9.6 KB

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

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

Download in other formats:

Original Format