source: LMDZ.3.3/trunk/libf/filtrez/filtreg.F @ 351

Last change on this file since 351 was 196, checked in by lmdz, 24 years ago

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