source: LMDZ4/branches/V3_test/libf/filtrez/filtreg.F @ 5297

Last change on this file since 5297 was 704, checked in by Laurent Fairhead, 18 years ago

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