Context Navigation

filtreg.F @ 5434

Last change on this file since 5434 was 1705, checked in by Ehouarn Millour, 12 years ago
Added arch files for ADA (IDRIS IBMx3750) and made the following code modifications: phylmd/printflag.F : removed "print" of unset variable (radpas0) dyn3dmem/integrd_loc.F : removed unecessary "include mpif.h" dyn3dmem/leapfrog_loc.F : removed unecessary "include mpif.h" and allocate saved variables at first call dyn3dmem/mod_filtreg_p.F : added matmul() alternatives to call to BLAS routine SGEMM (which was incorectly set as DGEMM; which would fail if running with -r4) filtrez/filtreg.F: changed calls to DGEMM into calls to SGEMM, so that code works with either -r4 or -r8 (the later being used in conjunction with "BLAS SGEMV=DGEMV SGEMM=DGEMM" preprocessing statements) EM
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 (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_sddv
115	sdd2_type = type_unsddv
116	ELSE
117	sdd1_type = type_unsddv
118	sdd2_type = type_sddv
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	#ifdef BLAS
187	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
188	& matrinvn(1,1,j),
189	& iim, champ(1,j,1), iip1*nlat, 0.0,
190	& eignq(1,j-jdfil+1,1), iim*nlat)
191	#else
192	eignq(:,j-jdfil+1,:)
193	$ = matmul(matrinvn(:,:,j), champ(:iim,j,:))
194	#endif
195	END DO
196
197	ELSE IF ( griscal ) THEN
198
199	DO j = jdfil,jffil
200	#ifdef BLAS
201	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
202	& matriceun(1,1,j),
203	& iim, champ(1,j,1), iip1*nlat, 0.0,
204	& eignq(1,j-jdfil+1,1), iim*nlat)
205	#else
206	eignq(:,j-jdfil+1,:)
207	$ = matmul(matriceun(:,:,j), champ(:iim,j,:))
208	#endif
209	END DO
210
211	ELSE
212
213	DO j = jdfil,jffil
214	#ifdef BLAS
215	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
216	& matricevn(1,1,j),
217	& iim, champ(1,j,1), iip1*nlat, 0.0,
218	& eignq(1,j-jdfil+1,1), iim*nlat)
219	#else
220	eignq(:,j-jdfil+1,:)
221	$ = matmul(matricevn(:,:,j), champ(:iim,j,:))
222	#endif
223	END DO
224
225	ENDIF
226
227	ELSE
228
229	IF( ifiltre. EQ. -2 ) THEN
230
231	DO j = jdfil,jffil
232	#ifdef BLAS
233	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
234	& matrinvs(1,1,j-jfiltsu+1),
235	& iim, champ(1,j,1), iip1*nlat, 0.0,
236	& eignq(1,j-jdfil+1,1), iim*nlat)
237	#else
238	eignq(:,j-jdfil+1,:)
239	$ = matmul(matrinvs(:,:,j-jfiltsu+1),
240	$ champ(:iim,j,:))
241	#endif
242	END DO
243
244
245	ELSE IF ( griscal ) THEN
246
247	DO j = jdfil,jffil
248	#ifdef BLAS
249	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
250	& matriceus(1,1,j-jfiltsu+1),
251	& iim, champ(1,j,1), iip1*nlat, 0.0,
252	& eignq(1,j-jdfil+1,1), iim*nlat)
253	#else
254	eignq(:,j-jdfil+1,:)
255	$ = matmul(matriceus(:,:,j-jfiltsu+1),
256	$ champ(:iim,j,:))
257	#endif
258	END DO
259
260	ELSE
261
262	DO j = jdfil,jffil
263	#ifdef BLAS
264	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0,
265	& matricevs(1,1,j-jfiltsv+1),
266	& iim, champ(1,j,1), iip1*nlat, 0.0,
267	& eignq(1,j-jdfil+1,1), iim*nlat)
268	#else
269	eignq(:,j-jdfil+1,:)
270	$ = matmul(matricevs(:,:,j-jfiltsv+1),
271	$ champ(:iim,j,:))
272	#endif
273	END DO
274
275	ENDIF
276
277	ENDIF
278
279	IF( ifiltre.EQ. 2 ) THEN
280
281	DO l = 1, nbniv
282	DO j = jdfil,jffil
283	DO i = 1, iim
284	champ( i,j,l ) =
285	& (champ(i,j,l) + eignq(i,j-jdfil+1,l))
286	& * sdd12(i,sdd2_type) ! sdd2(i)
287	END DO
288	END DO
289	END DO
290
291	ELSE
292
293	DO l = 1, nbniv
294	DO j = jdfil,jffil
295	DO i = 1, iim
296	champ( i,j,l ) =
297	& (champ(i,j,l) - eignq(i,j-jdfil+1,l))
298	& * sdd12(i,sdd2_type) ! sdd2(i)
299	END DO
300	END DO
301	END DO
302
303	ENDIF
304
305	DO l = 1, nbniv
306	DO j = jdfil,jffil
307	champ( iip1,j,l ) = champ( 1,j,l )
308	END DO
309	END DO
310
311
312	ENDDO
313
314	1111 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a
315	& filtrer, sur la grille des scalaires'/)
316	2222 FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
317	& ltrer, sur la grille de V ou de Z'/)
318	RETURN
319	END

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

Download in other formats:

Original Format