source: LMDZ.3.3/branches/rel-LF/libf/filtrez/filtreg.F @ 575

Last change on this file since 575 was 495, checked in by lmdzadmin, 21 years ago

Optimisation de differentes routines, IM, MAF, 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: 10.7 KB
Line 
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(2),jffil(2)
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)
58cIM   REAL  eignq(iim), sdd1(iim),sdd2(iim)
59      REAL  eignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim)
60      LOGICAL    griscal
61      INTEGER    hemisph, iaire
62c
63
64      IF(ifiltre.EQ.1.or.ifiltre.EQ.-1)
65     *    STOP'Pas de transformee simple dans cette version'
66
67      IF( iter.EQ. 2 )  THEN
68       PRINT *,' Pas d iteration du filtre dans cette version !'
69     * , ' Utiliser old_filtreg et repasser !'
70           STOP
71      ENDIF
72
73      IF( ifiltre.EQ. -2 .AND..NOT.griscal )     THEN
74       PRINT *,' Cette routine ne calcule le filtre inverse que ',
75     * ' sur la grille des scalaires !'
76           STOP
77      ENDIF
78
79      IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 )  THEN
80       PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
81     *,' corriger et repasser !'
82           STOP
83      ENDIF
84c
85
86      iim2   = iim * iim
87      immjm  = iim * jjm
88c
89c
90      IF( griscal )   THEN
91         IF( nlat. NE. jjp1 )  THEN
92             PRINT  1111
93             STOP
94         ELSE
95c
96             IF( iaire.EQ.1 )  THEN
97                CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 )
98                CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
99             ELSE
100                CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
101                CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
102             END IF
103c
104             jdfil1 = 2
105             jffil1 = jfiltnu
106             jdfil2 = jfiltsu
107             jffil2 = jjm
108          END IF
109      ELSE
110          IF( nlat.NE.jjm )  THEN
111             PRINT  2222
112             STOP
113          ELSE
114c
115             IF( iaire.EQ.1 )  THEN
116                CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 )
117                CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
118             ELSE
119                CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
120                CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
121             END IF
122c
123             jdfil1 = 1
124             jffil1 = jfiltnv
125             jdfil2 = jfiltsv
126             jffil2 = jjm
127          END IF
128      END IF
129c
130      jdfil(1) = jdfil1
131      jffil(1) = jffil1
132      jdfil(2) = jdfil2
133      jffil(2) = jffil2
134c
135      DO 100  hemisph = 1, 2
136c
137c     IF ( hemisph.EQ.1 )  THEN
138c         jdfil = jdfil1
139c         jffil = jffil1
140c     ELSE
141c         jdfil = jdfil2
142c         jffil = jffil2
143c     END IF
144
145 
146      DO 50  l = 1, nbniv
147      DO 30  j = jdfil(hemisph),jffil(hemisph)
148 
149 
150      DO  5  i = 1, iim
151      champ(i,j,l) = champ(i,j,l) * sdd1(i)
152   5  CONTINUE
153c
154 30   CONTINUE
155 50   CONTINUE
156
157      IF( hemisph. EQ. 1 )      THEN
158
159        IF( ifiltre. EQ. -2 )   THEN
160#ifdef CRAY
161      DO l = 1, nbniv
162      DO j = jdfil(hemisph),jffil(hemisph)
163         CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq  ,
164     *                             1, iim, iim                         )
165      ENDDO
166      ENDDO
167#else
168#ifdef BLAS
169      DO l = 1, nbniv
170      DO j = jdfil(hemisph),jffil(hemisph)
171      CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
172     .           champ(1,j,l), 1, 0.0, eignq, 1)
173      ENDDO
174      ENDDO
175#else
176      DO l = 1, nbniv
177      DO j = jdfil(hemisph),jffil(hemisph)
178      DO k = 1, iim
179c        eignq(k) = 0.0
180         eignq(k,j,l) = 0.0
181      ENDDO
182      DO k = 1, iim
183      DO i = 1, iim
184c        eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l)
185         eignq(k,j,l) = eignq(k,j,l) + matrinvn(k,i,j)*champ(i,j,l)
186      ENDDO
187      ENDDO
188      ENDDO
189      ENDDO
190#endif
191#endif
192        ELSE IF ( griscal )     THEN
193#ifdef CRAY
194      DO l = 1, nbniv
195      DO j = jdfil(hemisph),jffil(hemisph)
196         CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
197     *                             1, iim, iim                         )
198      ENDDO
199      ENDDO
200#else
201#ifdef BLAS
202      DO l = 1, nbniv
203      DO j = jdfil(hemisph),jffil(hemisph)
204      CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
205     .           champ(1,j,l), 1, 0.0, eignq, 1)
206      ENDDO
207      ENDDO
208#else
209      DO l = 1, nbniv
210      DO j = jdfil(hemisph),jffil(hemisph)
211      DO k = 1, iim
212c        eignq(k) = 0.0
213         eignq(k,j,l) = 0.0
214      ENDDO
215      DO i = 1, iim
216      DO k = 1, iim
217c        eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l)
218         eignq(k,j,l) = eignq(k,j,l) + matriceun(k,i,j)*champ(i,j,l)
219      ENDDO
220      ENDDO
221      ENDDO
222      ENDDO
223#endif
224#endif
225        ELSE
226#ifdef CRAY
227      DO l = 1, nbniv
228      DO j = jdfil(hemisph),jffil(hemisph)
229         CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
230     *                             1, iim, iim                         )
231      ENDDO
232      ENDDO
233#else
234#ifdef BLAS
235      DO l = 1, nbniv
236      DO j = jdfil(hemisph),jffil(hemisph)
237      CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
238     .           champ(1,j,l), 1, 0.0, eignq, 1)
239      ENDDO
240      ENDDO
241#else
242      DO l = 1, nbniv
243      DO j = jdfil(hemisph),jffil(hemisph)
244      DO k = 1, iim
245c        eignq(k) = 0.0
246         eignq(k,j,l) = 0.0
247      ENDDO
248      DO i = 1, iim
249      DO k = 1, iim
250c        eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l)
251         eignq(k,j,l) = eignq(k,j,l) + matricevn(k,i,j)*champ(i,j,l)
252      ENDDO
253      ENDDO
254      ENDDO
255      ENDDO
256#endif
257#endif
258        ENDIF
259
260      ELSE
261
262        IF( ifiltre. EQ. -2 )   THEN
263#ifdef CRAY
264      DO l = 1, nbniv
265      DO j = jdfil(hemisph),jffil(hemisph)
266         CALL MXVA( matrinvs(1,1,j-jfiltsu+1),  1, iim, champ(1,j,l),1 , 
267     *                          eignq,  1, iim, iim                    )
268      ENDDO
269      ENDDO
270#else
271#ifdef BLAS
272      DO l = 1, nbniv
273      DO j = jdfil(hemisph),jffil(hemisph)
274      CALL SGEMV("N", iim,iim, 1.0, matrinvs(1,1,j-jfiltsu+1),iim,
275     .           champ(1,j,l), 1, 0.0, eignq, 1)
276      ENDDO
277      ENDDO
278#else
279      DO l = 1, nbniv
280      DO j = jdfil(hemisph),jffil(hemisph)
281      DO k = 1, iim
282c        eignq(k) = 0.0
283         eignq(k,j,l) = 0.0
284      ENDDO
285      DO i = 1, iim
286      DO k = 1, iim
287c        eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l)
288         eignq(k,j,l) = eignq(k,j,l) +
289     .matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l)
290      ENDDO
291      ENDDO
292      ENDDO
293      ENDDO
294#endif
295#endif
296        ELSE IF ( griscal )     THEN
297#ifdef CRAY
298      DO l = 1, nbniv
299      DO j = jdfil(hemisph),jffil(hemisph)
300         CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
301     *                          eignq,  1, iim, iim                    )
302      ENDDO
303      ENDDO
304#else
305#ifdef BLAS
306      DO l = 1, nbniv
307      DO j = jdfil(hemisph),jffil(hemisph)
308      CALL SGEMV("N", iim,iim, 1.0, matriceus(1,1,j-jfiltsu+1),iim,
309     .           champ(1,j,l), 1, 0.0, eignq, 1)
310      ENDDO
311      ENDDO
312#else
313      DO l = 1, nbniv
314      DO j = jdfil(hemisph),jffil(hemisph)
315      DO k = 1, iim
316c        eignq(k) = 0.0
317         eignq(k,j,l) = 0.0
318      ENDDO
319      DO i = 1, iim
320      DO k = 1, iim
321c        eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l)
322         eignq(k,j,l) = eignq(k,j,l) +
323     .matriceus(k,i,j-jfiltsu+1)*champ(i,j,l)
324      ENDDO
325      ENDDO
326      ENDDO
327      ENDDO
328#endif
329#endif
330        ELSE
331#ifdef CRAY
332      DO l = 1, nbniv
333      DO j = jdfil(hemisph),jffil(hemisph)
334         CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 ,
335     *                          eignq,  1, iim, iim                    )
336      ENDDO
337      ENDDO
338#else
339#ifdef BLAS
340      DO l = 1, nbniv
341      DO j = jdfil(hemisph),jffil(hemisph)
342      CALL SGEMV("N", iim,iim, 1.0, matricevs(1,1,j-jfiltsv+1),iim,
343     .           champ(1,j,l), 1, 0.0, eignq, 1)
344      ENDDO
345      ENDDO
346#else
347      DO l = 1, nbniv
348      DO j = jdfil(hemisph),jffil(hemisph)
349      DO k = 1, iim
350c        eignq(k) = 0.0
351         eignq(k,j,l) = 0.0
352      ENDDO
353      DO i = 1, iim
354      DO k = 1, iim
355c        eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l)
356         eignq(k,j,l) = eignq(k,j,l) +
357     .matricevs(k,i,j-jfiltsv+1)*champ(i,j,l)
358      ENDDO
359      ENDDO
360      ENDDO
361      ENDDO
362#endif
363#endif
364        ENDIF
365
366      ENDIF
367c
368      IF( ifiltre.EQ. 2 )  THEN
369      DO l = 1, nbniv
370      DO j = jdfil(hemisph),jffil(hemisph)
371        DO 15 i = 1, iim
372c       champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
373        champ( i,j,l ) = ( champ(i,j,l) + eignq(i,j,l) ) * sdd2(i)
374  15    CONTINUE
375      ENDDO
376      ENDDO
377      ELSE
378      DO l = 1, nbniv
379      DO j = jdfil(hemisph),jffil(hemisph)
380        DO 16 i=1,iim
381c       champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
382        champ( i,j,l ) = ( champ(i,j,l) - eignq(i,j,l) ) * sdd2(i)
38316      CONTINUE
384      ENDDO
385      ENDDO
386      ENDIF
387c
388      DO l = 1, nbniv
389      DO j = jdfil(hemisph),jffil(hemisph)
390      champ( iip1,j,l ) = champ( 1,j,l )
391      ENDDO
392      ENDDO
393c
394c 30  CONTINUE
395c
396c 50  CONTINUE
397c   
398 100  CONTINUE
399c
4001111  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a
401     *filtrer, sur la grille des scalaires'/)
4022222  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
403     *ltrer, sur la grille de V ou de Z'/)
404      RETURN
405      END
Note: See TracBrowser for help on using the repository browser.