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
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
[495]51      INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil(2),jffil(2)
[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)
[495]58cIM   REAL  eignq(iim), sdd1(iim),sdd2(iim)
59      REAL  eignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim)
[2]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
[495]130      jdfil(1) = jdfil1
131      jffil(1) = jffil1
132      jdfil(2) = jdfil2
133      jffil(2) = jffil2
[2]134c
135      DO 100  hemisph = 1, 2
136c
[495]137c     IF ( hemisph.EQ.1 )  THEN
138c         jdfil = jdfil1
139c         jffil = jffil1
140c     ELSE
141c         jdfil = jdfil2
142c         jffil = jffil2
143c     END IF
[2]144
145 
146      DO 50  l = 1, nbniv
[495]147      DO 30  j = jdfil(hemisph),jffil(hemisph)
[2]148 
149 
150      DO  5  i = 1, iim
151      champ(i,j,l) = champ(i,j,l) * sdd1(i)
152   5  CONTINUE
153c
[495]154 30   CONTINUE
155 50   CONTINUE
[2]156
157      IF( hemisph. EQ. 1 )      THEN
158
159        IF( ifiltre. EQ. -2 )   THEN
160#ifdef CRAY
[495]161      DO l = 1, nbniv
162      DO j = jdfil(hemisph),jffil(hemisph)
[2]163         CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq  ,
164     *                             1, iim, iim                         )
[495]165      ENDDO
166      ENDDO
[2]167#else
168#ifdef BLAS
[495]169      DO l = 1, nbniv
170      DO j = jdfil(hemisph),jffil(hemisph)
[2]171      CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
172     .           champ(1,j,l), 1, 0.0, eignq, 1)
[495]173      ENDDO
174      ENDDO
[2]175#else
[495]176      DO l = 1, nbniv
177      DO j = jdfil(hemisph),jffil(hemisph)
[2]178      DO k = 1, iim
[495]179c        eignq(k) = 0.0
180         eignq(k,j,l) = 0.0
[2]181      ENDDO
182      DO k = 1, iim
183      DO i = 1, iim
[495]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)
[2]186      ENDDO
187      ENDDO
[495]188      ENDDO
189      ENDDO
[2]190#endif
191#endif
192        ELSE IF ( griscal )     THEN
193#ifdef CRAY
[495]194      DO l = 1, nbniv
195      DO j = jdfil(hemisph),jffil(hemisph)
[2]196         CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
197     *                             1, iim, iim                         )
[495]198      ENDDO
199      ENDDO
[2]200#else
201#ifdef BLAS
[495]202      DO l = 1, nbniv
203      DO j = jdfil(hemisph),jffil(hemisph)
[2]204      CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
205     .           champ(1,j,l), 1, 0.0, eignq, 1)
[495]206      ENDDO
207      ENDDO
[2]208#else
[495]209      DO l = 1, nbniv
210      DO j = jdfil(hemisph),jffil(hemisph)
[2]211      DO k = 1, iim
[495]212c        eignq(k) = 0.0
213         eignq(k,j,l) = 0.0
[2]214      ENDDO
215      DO i = 1, iim
216      DO k = 1, iim
[495]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)
[2]219      ENDDO
220      ENDDO
[495]221      ENDDO
222      ENDDO
[2]223#endif
224#endif
225        ELSE
226#ifdef CRAY
[495]227      DO l = 1, nbniv
228      DO j = jdfil(hemisph),jffil(hemisph)
[2]229         CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
230     *                             1, iim, iim                         )
[495]231      ENDDO
232      ENDDO
[2]233#else
234#ifdef BLAS
[495]235      DO l = 1, nbniv
236      DO j = jdfil(hemisph),jffil(hemisph)
[2]237      CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
238     .           champ(1,j,l), 1, 0.0, eignq, 1)
[495]239      ENDDO
240      ENDDO
[2]241#else
[495]242      DO l = 1, nbniv
243      DO j = jdfil(hemisph),jffil(hemisph)
[2]244      DO k = 1, iim
[495]245c        eignq(k) = 0.0
246         eignq(k,j,l) = 0.0
[2]247      ENDDO
248      DO i = 1, iim
249      DO k = 1, iim
[495]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)
[2]252      ENDDO
253      ENDDO
[495]254      ENDDO
255      ENDDO
[2]256#endif
257#endif
258        ENDIF
259
260      ELSE
261
262        IF( ifiltre. EQ. -2 )   THEN
263#ifdef CRAY
[495]264      DO l = 1, nbniv
265      DO j = jdfil(hemisph),jffil(hemisph)
[2]266         CALL MXVA( matrinvs(1,1,j-jfiltsu+1),  1, iim, champ(1,j,l),1 , 
267     *                          eignq,  1, iim, iim                    )
[495]268      ENDDO
269      ENDDO
[2]270#else
271#ifdef BLAS
[495]272      DO l = 1, nbniv
273      DO j = jdfil(hemisph),jffil(hemisph)
[2]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)
[495]276      ENDDO
277      ENDDO
[2]278#else
[495]279      DO l = 1, nbniv
280      DO j = jdfil(hemisph),jffil(hemisph)
[2]281      DO k = 1, iim
[495]282c        eignq(k) = 0.0
283         eignq(k,j,l) = 0.0
[2]284      ENDDO
285      DO i = 1, iim
286      DO k = 1, iim
[495]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)
[2]290      ENDDO
291      ENDDO
[495]292      ENDDO
293      ENDDO
[2]294#endif
295#endif
296        ELSE IF ( griscal )     THEN
297#ifdef CRAY
[495]298      DO l = 1, nbniv
299      DO j = jdfil(hemisph),jffil(hemisph)
[2]300         CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 ,
301     *                          eignq,  1, iim, iim                    )
[495]302      ENDDO
303      ENDDO
[2]304#else
305#ifdef BLAS
[495]306      DO l = 1, nbniv
307      DO j = jdfil(hemisph),jffil(hemisph)
[2]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)
[495]310      ENDDO
311      ENDDO
[2]312#else
[495]313      DO l = 1, nbniv
314      DO j = jdfil(hemisph),jffil(hemisph)
[2]315      DO k = 1, iim
[495]316c        eignq(k) = 0.0
317         eignq(k,j,l) = 0.0
[2]318      ENDDO
319      DO i = 1, iim
320      DO k = 1, iim
[495]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)
[2]324      ENDDO
325      ENDDO
[495]326      ENDDO
327      ENDDO
[2]328#endif
329#endif
330        ELSE
331#ifdef CRAY
[495]332      DO l = 1, nbniv
333      DO j = jdfil(hemisph),jffil(hemisph)
[2]334         CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 ,
335     *                          eignq,  1, iim, iim                    )
[495]336      ENDDO
337      ENDDO
[2]338#else
339#ifdef BLAS
[495]340      DO l = 1, nbniv
341      DO j = jdfil(hemisph),jffil(hemisph)
[2]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)
[495]344      ENDDO
345      ENDDO
[2]346#else
[495]347      DO l = 1, nbniv
348      DO j = jdfil(hemisph),jffil(hemisph)
[2]349      DO k = 1, iim
[495]350c        eignq(k) = 0.0
351         eignq(k,j,l) = 0.0
[2]352      ENDDO
353      DO i = 1, iim
354      DO k = 1, iim
[495]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)
[2]358      ENDDO
359      ENDDO
[495]360      ENDDO
361      ENDDO
[2]362#endif
363#endif
364        ENDIF
365
366      ENDIF
367c
368      IF( ifiltre.EQ. 2 )  THEN
[495]369      DO l = 1, nbniv
370      DO j = jdfil(hemisph),jffil(hemisph)
[2]371        DO 15 i = 1, iim
[495]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)
[2]374  15    CONTINUE
[495]375      ENDDO
376      ENDDO
[2]377      ELSE
[495]378      DO l = 1, nbniv
379      DO j = jdfil(hemisph),jffil(hemisph)
[2]380        DO 16 i=1,iim
[495]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)
[2]38316      CONTINUE
[495]384      ENDDO
385      ENDDO
[2]386      ENDIF
387c
[495]388      DO l = 1, nbniv
389      DO j = jdfil(hemisph),jffil(hemisph)
[2]390      champ( iip1,j,l ) = champ( 1,j,l )
[495]391      ENDDO
392      ENDDO
[2]393c
[495]394c 30  CONTINUE
[2]395c
[495]396c 50  CONTINUE
[2]397c   
398 100  CONTINUE
399c
[231]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'/)
[2]404      RETURN
405      END
Note: See TracBrowser for help on using the repository browser.