Index: LMDZ4/branches/LMDZ4-dev/libf/filtrez/filtreg.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/filtrez/filtreg.F	(revision 1074)
+++ LMDZ4/branches/LMDZ4-dev/libf/filtrez/filtreg.F	(revision 1086)
@@ -3,6 +3,8 @@
 !
       SUBROUTINE filtreg ( champ, nlat, nbniv, ifiltre,iaire,
-     .   griscal ,iter)
-
+     &     griscal ,iter)
+      
+      USE filtre
+      
       IMPLICIT NONE
 c=======================================================================
@@ -46,255 +48,239 @@
 #include "dimensions.h"
 #include "paramet.h"
-#include "parafilt.h"
 #include "coefils.h"
-c
-      INTEGER nlat,nbniv,ifiltre,iter
-      INTEGER i,j,l,k
-      INTEGER iim2,immjm
-      INTEGER jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
-
-      REAL  champ( iip1,nlat,nbniv)
-      REAL matriceun,matriceus,matricevn,matricevs,matrinvn,matrinvs
-      COMMON/matrfil/matriceun(iim,iim,nfilun),matriceus(iim,iim,nfilus)
-     ,             , matricevn(iim,iim,nfilvn),matricevs(iim,iim,nfilvs)
-     ,             ,  matrinvn(iim,iim,nfilun),matrinvs (iim,iim,nfilus)
-      REAL  eignq(iim), sdd1(iim),sdd2(iim)
+
+      INTEGER    nlat,nbniv,ifiltre,iter
+      INTEGER    i,j,l,k
+      INTEGER    iim2,immjm
+      INTEGER    jdfil1,jdfil2,jffil1,jffil2,jdfil,jffil
+
+      REAL       champ( iip1,nlat,nbniv)
+
+      REAL       eignq(iim,nlat,nbniv), sdd1(iim),sdd2(iim)
       LOGICAL    griscal
       INTEGER    hemisph, iaire
-c
+
+      LOGICAL,SAVE     :: first=.TRUE.
+
+      REAL, SAVE :: sdd12(iim,4)
+
+      INTEGER, PARAMETER :: type_sddu=1
+      INTEGER, PARAMETER :: type_sddv=2
+      INTEGER, PARAMETER :: type_unsddu=3
+      INTEGER, PARAMETER :: type_unsddv=4
+
+      INTEGER :: sdd1_type, sdd2_type
+
+      IF (first) THEN
+         sdd12(1:iim,type_sddu) = sddu(1:iim)
+         sdd12(1:iim,type_sddv) = sddv(1:iim)
+         sdd12(1:iim,type_unsddu) = unsddu(1:iim)
+         sdd12(1:iim,type_unsddv) = unsddv(1:iim)
+
+         first=.FALSE.
+      ENDIF
 
       IF(ifiltre.EQ.1.or.ifiltre.EQ.-1) 
-     *    STOP'Pas de transformee simple dans cette version'
-
+     &     STOP'Pas de transformee simple dans cette version'
+      
       IF( iter.EQ. 2 )  THEN
-       PRINT *,' Pas d iteration du filtre dans cette version !'
-     * , ' Utiliser old_filtreg et repasser !'
-           STOP
-      ENDIF
-
+         PRINT *,' Pas d iteration du filtre dans cette version !'
+     &        , ' Utiliser old_filtreg et repasser !'
+         STOP
+      ENDIF
+      
       IF( ifiltre.EQ. -2 .AND..NOT.griscal )     THEN
-       PRINT *,' Cette routine ne calcule le filtre inverse que ',
-     * ' sur la grille des scalaires !'
-           STOP
-      ENDIF
-
+         PRINT *,' Cette routine ne calcule le filtre inverse que '
+     &        , ' sur la grille des scalaires !'
+         STOP
+      ENDIF
+      
       IF( ifiltre.NE.2 .AND.ifiltre.NE. - 2 )  THEN
-       PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
-     *,' corriger et repasser !'
-           STOP
-      ENDIF
-c
-
+         PRINT *,' Probleme dans filtreg car ifiltre NE 2 et NE -2'
+     &        , ' corriger et repasser !'
+         STOP
+      ENDIF
+      
       iim2   = iim * iim
       immjm  = iim * jjm
-c
-c
+
       IF( griscal )   THEN
          IF( nlat. NE. jjp1 )  THEN
-             PRINT  1111
-             STOP
-         ELSE
-c
-             IF( iaire.EQ.1 )  THEN
-                CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 ) 
-                CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
-             ELSE
-                CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
-                CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
-             END IF
-c
-             jdfil1 = 2
-             jffil1 = jfiltnu
-             jdfil2 = jfiltsu
-             jffil2 = jjm
-          END IF
+            PRINT  1111
+            STOP
+         ELSE
+            
+            IF( iaire.EQ.1 )  THEN
+               sdd1_type = type_sddu
+               sdd2_type = type_unsddu
+            ELSE
+               sdd1_type = type_unsddu
+               sdd2_type = type_sddu
+            ENDIF
+
+c            IF( iaire.EQ.1 )  THEN
+c               CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 ) 
+c               CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
+c            ELSE
+c               CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
+c               CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
+c            END IF
+            
+            jdfil1 = 2
+            jffil1 = jfiltnu
+            jdfil2 = jfiltsu
+            jffil2 = jjm
+         END IF
       ELSE
-          IF( nlat.NE.jjm )  THEN
-             PRINT  2222
-             STOP
-          ELSE
-c
-             IF( iaire.EQ.1 )  THEN
-                CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 ) 
-                CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
-             ELSE
-                CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
-                CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
-             END IF
-c
-             jdfil1 = 1
-             jffil1 = jfiltnv
-             jdfil2 = jfiltsv
-             jffil2 = jjm
-          END IF
+         IF( nlat.NE.jjm )  THEN
+            PRINT  2222
+            STOP
+         ELSE
+            
+            IF( iaire.EQ.1 )  THEN
+               sdd1_type = type_sddu
+               sdd2_type = type_unsddu
+            ELSE
+               sdd1_type = type_unsddu
+               sdd2_type = type_sddu
+            ENDIF
+
+c            IF( iaire.EQ.1 )  THEN
+c               CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 ) 
+c               CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
+c            ELSE
+c               CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
+c               CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
+c            END IF
+            
+            jdfil1 = 1
+            jffil1 = jfiltnv
+            jdfil2 = jfiltsv
+            jffil2 = jjm
+         END IF
       END IF
-c
-c
-      DO 100  hemisph = 1, 2
-c
-      IF ( hemisph.EQ.1 )  THEN
-          jdfil = jdfil1
-          jffil = jffil1
-      ELSE
-          jdfil = jdfil2
-          jffil = jffil2
-      END IF
-
- 
-      DO 50  l = 1, nbniv
-      DO 30  j = jdfil,jffil
- 
- 
-      DO  5  i = 1, iim
-      champ(i,j,l) = champ(i,j,l) * sdd1(i)
-   5  CONTINUE
-c
-
-      IF( hemisph. EQ. 1 )      THEN
-
-        IF( ifiltre. EQ. -2 )   THEN
-#ifdef CRAY
-         CALL MXVA( matrinvn(1,1,j), 1, iim, champ(1,j,l), 1, eignq  , 
-     *                             1, iim, iim                         )
-#else
-#ifdef BLAS
-      CALL SGEMV("N", iim,iim, 1.0, matrinvn(1,1,j),iim,
-     .           champ(1,j,l), 1, 0.0, eignq, 1)
-#else
-      DO k = 1, iim
-         eignq(k) = 0.0
+      
+      DO hemisph = 1, 2
+         
+         IF ( hemisph.EQ.1 )  THEN
+            jdfil = jdfil1
+            jffil = jffil1
+         ELSE
+            jdfil = jdfil2
+            jffil = jffil2
+         END IF
+         
+         DO l = 1, nbniv
+            DO j = jdfil,jffil
+               DO i = 1, iim
+                  champ(i,j,l) = champ(i,j,l) * sdd12(i,sdd1_type) ! sdd1(i)
+               END DO
+            END DO
+         END DO
+         
+         IF( hemisph. EQ. 1 )      THEN
+            
+            IF( ifiltre. EQ. -2 )   THEN
+               
+               DO j = jdfil,jffil
+                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 
+     &                 matrinvn(1,1,j),
+     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
+     &                 eignq(1,j-jdfil+1,1), iim*nlat)
+               END DO
+               
+            ELSE IF ( griscal )     THEN
+               
+               DO j = jdfil,jffil
+                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 
+     &                 matriceun(1,1,j),
+     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
+     &                 eignq(1,j-jdfil+1,1), iim*nlat)
+               END DO
+               
+            ELSE 
+               
+               DO j = jdfil,jffil
+                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 
+     &                 matricevn(1,1,j),
+     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
+     &                 eignq(1,j-jdfil+1,1), iim*nlat)
+               END DO
+               
+            ENDIF
+            
+         ELSE
+            
+            IF( ifiltre. EQ. -2 )   THEN
+               
+               DO j = jdfil,jffil
+                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 
+     &                 matrinvs(1,1,j-jfiltsu+1),
+     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
+     &                 eignq(1,j-jdfil+1,1), iim*nlat)
+               END DO
+               
+               
+            ELSE IF ( griscal )     THEN
+               
+               DO j = jdfil,jffil
+                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 
+     &                 matriceus(1,1,j-jfiltsu+1),
+     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
+     &                 eignq(1,j-jdfil+1,1), iim*nlat)
+               END DO
+                              
+            ELSE 
+               
+               DO j = jdfil,jffil
+                  CALL DGEMM("N", "N", iim, nbniv, iim, 1.0, 
+     &                 matricevs(1,1,j-jfiltsv+1),
+     &                 iim, champ(1,j,1), iip1*nlat, 0.0,
+     &                 eignq(1,j-jdfil+1,1), iim*nlat)
+               END DO
+                              
+            ENDIF
+            
+         ENDIF
+         
+         IF( ifiltre.EQ. 2 )  THEN
+            
+            DO l = 1, nbniv
+               DO j = jdfil,jffil
+                  DO i = 1, iim
+                     champ( i,j,l ) = 
+     &                    (champ(i,j,l) + eignq(i,j-jdfil+1,l))
+     &                    * sdd12(i,sdd2_type) ! sdd2(i)
+                  END DO
+               END DO
+            END DO
+
+         ELSE
+
+            DO l = 1, nbniv
+               DO j = jdfil,jffil
+                  DO i = 1, iim
+                     champ( i,j,l ) = 
+     &                    (champ(i,j,l) - eignq(i,j-jdfil+1,l))
+     &                    * sdd12(i,sdd2_type) ! sdd2(i)
+                  END DO
+               END DO
+            END DO
+
+         ENDIF
+
+         DO l = 1, nbniv
+            DO j = jdfil,jffil
+               champ( iip1,j,l ) = champ( 1,j,l )
+            END DO
+         END DO
+
+     
       ENDDO
-      DO k = 1, iim
-      DO i = 1, iim
-         eignq(k) = eignq(k) + matrinvn(k,i,j)*champ(i,j,l)
-      ENDDO
-      ENDDO
-#endif
-#endif
-        ELSE IF ( griscal )     THEN
-#ifdef CRAY
-         CALL MXVA( matriceun(1,1,j), 1, iim, champ(1,j,l), 1, eignq ,
-     *                             1, iim, iim                         )
-#else
-#ifdef BLAS
-      CALL SGEMV("N", iim,iim, 1.0, matriceun(1,1,j),iim,
-     .           champ(1,j,l), 1, 0.0, eignq, 1)
-#else
-      DO k = 1, iim
-         eignq(k) = 0.0
-      ENDDO
-      DO i = 1, iim
-      DO k = 1, iim
-         eignq(k) = eignq(k) + matriceun(k,i,j)*champ(i,j,l)
-      ENDDO
-      ENDDO
-#endif
-#endif
-        ELSE 
-#ifdef CRAY
-         CALL MXVA( matricevn(1,1,j), 1, iim, champ(1,j,l), 1, eignq , 
-     *                             1, iim, iim                         )
-#else
-#ifdef BLAS
-      CALL SGEMV("N", iim,iim, 1.0, matricevn(1,1,j),iim,
-     .           champ(1,j,l), 1, 0.0, eignq, 1)
-#else
-      DO k = 1, iim
-         eignq(k) = 0.0
-      ENDDO
-      DO i = 1, iim
-      DO k = 1, iim
-         eignq(k) = eignq(k) + matricevn(k,i,j)*champ(i,j,l)
-      ENDDO
-      ENDDO
-#endif
-#endif
-        ENDIF
-
-      ELSE
-
-        IF( ifiltre. EQ. -2 )   THEN
-#ifdef CRAY
-         CALL MXVA( matrinvs(1,1,j-jfiltsu+1),  1, iim, champ(1,j,l),1 ,  
-     *                          eignq,  1, iim, iim                    )
-#else
-#ifdef BLAS
-      CALL SGEMV("N", iim,iim, 1.0, matrinvs(1,1,j-jfiltsu+1),iim,
-     .           champ(1,j,l), 1, 0.0, eignq, 1)
-#else
-      DO k = 1, iim
-         eignq(k) = 0.0
-      ENDDO
-      DO i = 1, iim
-      DO k = 1, iim
-         eignq(k) = eignq(k) + matrinvs(k,i,j-jfiltsu+1)*champ(i,j,l)
-      ENDDO
-      ENDDO
-#endif
-#endif
-        ELSE IF ( griscal )     THEN
-#ifdef CRAY
-         CALL MXVA( matriceus(1,1,j-jfiltsu+1), 1, iim, champ(1,j,l),1 , 
-     *                          eignq,  1, iim, iim                    )
-#else
-#ifdef BLAS
-      CALL SGEMV("N", iim,iim, 1.0, matriceus(1,1,j-jfiltsu+1),iim,
-     .           champ(1,j,l), 1, 0.0, eignq, 1)
-#else
-      DO k = 1, iim
-         eignq(k) = 0.0
-      ENDDO
-      DO i = 1, iim
-      DO k = 1, iim
-         eignq(k) = eignq(k) + matriceus(k,i,j-jfiltsu+1)*champ(i,j,l)
-      ENDDO
-      ENDDO
-#endif
-#endif
-        ELSE 
-#ifdef CRAY
-         CALL MXVA( matricevs(1,1,j-jfiltsv+1), 1, iim, champ(1,j,l),1 , 
-     *                          eignq,  1, iim, iim                    )
-#else
-#ifdef BLAS
-      CALL SGEMV("N", iim,iim, 1.0, matricevs(1,1,j-jfiltsv+1),iim,
-     .           champ(1,j,l), 1, 0.0, eignq, 1)
-#else
-      DO k = 1, iim
-         eignq(k) = 0.0
-      ENDDO
-      DO i = 1, iim
-      DO k = 1, iim
-         eignq(k) = eignq(k) + matricevs(k,i,j-jfiltsv+1)*champ(i,j,l)
-      ENDDO
-      ENDDO
-#endif
-#endif
-        ENDIF
-
-      ENDIF
-c
-      IF( ifiltre.EQ. 2 )  THEN
-        DO 15 i = 1, iim
-        champ( i,j,l ) = ( champ(i,j,l) + eignq(i) ) * sdd2(i)
-  15    CONTINUE
-      ELSE
-        DO 16 i=1,iim
-        champ( i,j,l ) = ( champ(i,j,l) - eignq(i) ) * sdd2(i)
-16      CONTINUE
-      ENDIF
-c
-      champ( iip1,j,l ) = champ( 1,j,l )
-c
-  30  CONTINUE
-c
-  50  CONTINUE
-c    
- 100  CONTINUE
-c
+
 1111  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau  CHAMP a 
-     *filtrer, sur la grille des scalaires'/)
+     &     filtrer, sur la grille des scalaires'/)
 2222  FORMAT(//20x,'ERREUR dans le dimensionnement du tableau CHAMP a fi
-     *ltrer, sur la grille de V ou de Z'/)
+     &     ltrer, sur la grille de V ou de Z'/)
       RETURN
       END
