My Project
 All Classes Files Functions Variables Macros
filtreg.F
Go to the documentation of this file.
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,
315 'ERREUR dans le dimensionnement du tableau CHAMP a & filtrer, sur la grille des scalaires'/)
316 2222 FORMAT(//20x,
317 'ERREUR dans le dimensionnement du tableau CHAMP a fi & ltrer, sur la grille de V ou de Z'/)
318  RETURN
319  END