Index: LMDZ4/branches/LMDZ4-dev/libf/filtrez/coefils.h
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/filtrez/coefils.h	(revision 1074)
+++ LMDZ4/branches/LMDZ4-dev/libf/filtrez/coefils.h	(revision 1086)
@@ -2,9 +2,9 @@
 ! $Header$
 !
-      COMMON/coefils/jfiltnu,jfiltsu,jfiltnv,jfiltsv,sddu(iim),sddv(iim)
-     * ,unsddu(iim),unsddv(iim),coefilu(iim,jjm),coefilv(iim,jjm),
-     * modfrstu(jjm),modfrstv(jjm),eignfnu(iim,iim),eignfnv(iim,iim)
-     * ,coefilu2(iim,jjm),coefilv2(iim,jjm)
-c
+      COMMON/coefils/jfiltnu,jfiltsu,jfiltnv,jfiltsv,sddu(iim),sddv(iim)&
+     & ,unsddu(iim),unsddv(iim),coefilu(iim,jjm),coefilv(iim,jjm),      &
+     & modfrstu(jjm),modfrstv(jjm),eignfnu(iim,iim),eignfnv(iim,iim)    &
+     & ,coefilu2(iim,jjm),coefilv2(iim,jjm)
+!c
       INTEGER jfiltnu,jfiltsu,jfiltnv,jfiltsv,modfrstu,modfrstv
       REAL    sddu,sddv,unsddu,unsddv,coefilu,coefilv,eignfnu,eignfnv
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
Index: LMDZ4/branches/LMDZ4-dev/libf/filtrez/inifgn.F
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/filtrez/inifgn.F	(revision 1074)
+++ LMDZ4/branches/LMDZ4-dev/libf/filtrez/inifgn.F	(revision 1086)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Header: /home/cvsroot/LMDZ4/libf/filtrez/inifgn.F,v 1.1.1.1 2004-05-19 12:53:09 lmdzadmin Exp $
 !
       SUBROUTINE inifgn(dv)
Index: LMDZ4/branches/LMDZ4-dev/libf/filtrez/module_filtreg.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/filtrez/module_filtreg.F90	(revision 1086)
+++ LMDZ4/branches/LMDZ4-dev/libf/filtrez/module_filtreg.F90	(revision 1086)
@@ -0,0 +1,530 @@
+MODULE filtre
+
+  REAL, DIMENSION(:,:,:), ALLOCATABLE :: matriceun,matriceus,matricevn
+  REAL, DIMENSION(:,:,:), ALLOCATABLE :: matricevs,matrinvn,matrinvs
+
+CONTAINS
+
+  SUBROUTINE inifilr
+    !
+    !    ... H. Upadhyaya, O.Sharma   ...
+    !
+    IMPLICIT NONE
+    !
+    !     version 3 .....
+
+    !     Correction  le 28/10/97    P. Le Van .
+    !  -------------------------------------------------------------------
+#include "dimensions.h"
+#include "paramet.h"
+    !  -------------------------------------------------------------------
+#include "comgeom.h"
+#include "coefils.h"
+#include "logic.h"
+#include "serre.h"
+
+    REAL  dlonu(iim),dlatu(jjm)
+    REAL  rlamda( iim ),  eignvl( iim )
+    !
+
+    REAL    lamdamax,pi,cof
+    INTEGER i,j,modemax,imx,k,kf,ii
+    REAL dymin,dxmin,colat0
+    REAL eignft(iim,iim), coff
+
+    LOGICAL, SAVE :: first_call_inifilr = .TRUE.
+
+#ifdef CRAY
+    INTEGER   ISMIN
+    EXTERNAL  ISMIN
+    INTEGER iymin 
+    INTEGER ixmineq
+#endif
+    EXTERNAL  inifgn
+    !
+    ! ------------------------------------------------------------
+    !   This routine computes the eigenfunctions of the laplacien
+    !   on the stretched grid, and the filtering coefficients
+    !      
+    !  We designate:
+    !   eignfn   eigenfunctions of the discrete laplacien
+    !   eigenvl  eigenvalues
+    !   jfiltn   indexof the last scalar line filtered in NH
+    !   jfilts   index of the first line filtered in SH
+    !   modfrst  index of the mode from WHERE modes are filtered
+    !   modemax  maximum number of modes ( im )
+    !   coefil   filtering coefficients ( lamda_max*COS(rlat)/lamda )
+    !   sdd      SQRT( dx )
+    !      
+    !     the modes are filtered from modfrst to modemax
+    !      
+    !-----------------------------------------------------------
+    !
+
+    pi       = 2. * ASIN( 1. )
+
+    DO i = 1,iim
+       dlonu(i) = xprimu( i )
+    ENDDO
+    !
+    CALL inifgn(eignvl)
+    !
+    PRINT *,' EIGNVL '
+    PRINT 250,eignvl
+250 FORMAT( 1x,5e13.6)
+    !
+    ! compute eigenvalues and eigenfunctions
+    !
+    !
+    !.................................................................
+    !
+    !  compute the filtering coefficients for scalar lines and 
+    !  meridional wind v-lines
+    !
+    !  we filter all those latitude lines WHERE coefil < 1
+    !  NO FILTERING AT POLES
+    !
+    !  colat0 is to be used  when alpha (stretching coefficient)
+    !  is set equal to zero for the regular grid CASE 
+    !
+    !    .......   Calcul  de  colat0   .........
+    !     .....  colat0 = minimum de ( 0.5, min dy/ min dx )   ...
+    !
+    !
+    DO j = 1,jjm
+       dlatu( j ) = rlatu( j ) - rlatu( j+1 )
+    ENDDO
+    !
+#ifdef CRAY
+    iymin   = ISMIN( jjm, dlatu, 1 )
+    ixmineq = ISMIN( iim, dlonu, 1 )
+    dymin   = dlatu( iymin )
+    dxmin   = dlonu( ixmineq )
+#else
+    dxmin   =  dlonu(1)
+    DO  i  = 2, iim
+       dxmin = MIN( dxmin,dlonu(i) )
+    ENDDO
+    dymin  = dlatu(1)
+    DO j  = 2, jjm
+       dymin = MIN( dymin,dlatu(j) )
+    ENDDO
+#endif
+    !
+    !
+    colat0  =  MIN( 0.5, dymin/dxmin )
+    !
+    IF( .NOT.fxyhypb.AND.ysinus )  THEN
+       colat0 = 0.6
+       !         ...... a revoir  pour  ysinus !   .......
+       alphax = 0.
+    ENDIF
+    !
+    PRINT 50, colat0,alphax
+50  FORMAT(/15x,' Inifilr colat0 alphax ',2e16.7)
+    !
+    IF(alphax.EQ.1. )  THEN
+       PRINT *,' Inifilr  alphax doit etre  <  a 1.  Corriger '
+       STOP
+    ENDIF
+    !
+    lamdamax = iim / ( pi * colat0 * ( 1. - alphax ) )
+
+    !                        ... Correction  le 28/10/97  ( P.Le Van ) ..
+    !
+    DO i = 2,iim
+       rlamda( i ) = lamdamax/ SQRT( ABS( eignvl(i) ) )
+    ENDDO
+    !
+
+    DO j = 1,jjm
+       DO i = 1,iim
+          coefilu( i,j )  = 0.0
+          coefilv( i,j )  = 0.0
+          coefilu2( i,j ) = 0.0
+          coefilv2( i,j ) = 0.0
+       ENDDO
+    ENDDO
+
+    !
+    !    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....
+    !    .........................................................
+    !
+    modemax = iim
+
+!!!!    imx = modemax - 4 * (modemax/iim)
+
+    imx  = iim
+    !
+    PRINT *,' TRUNCATION AT ',imx
+    !
+    DO j = 2, jjm/2+1
+       cof = COS( rlatu(j) )/ colat0
+       IF ( cof .LT. 1. ) THEN
+          IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) jfiltnu= j
+       ENDIF
+
+       cof = COS( rlatu(jjp1-j+1) )/ colat0
+       IF ( cof .LT. 1. ) THEN
+          IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) &
+               jfiltsu= jjp1-j+1
+       ENDIF
+    ENDDO
+    !
+    DO j = 1, jjm/2
+       cof = COS( rlatv(j) )/ colat0
+       IF ( cof .LT. 1. ) THEN
+          IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) jfiltnv= j
+       ENDIF
+
+       cof = COS( rlatv(jjm-j+1) )/ colat0
+       IF ( cof .LT. 1. ) THEN
+          IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) &
+               jfiltsv= jjm-j+1
+       ENDIF
+    ENDDO
+    !                                 
+
+    IF ( jfiltnu.LE.0 ) jfiltnu=1
+    IF( jfiltnu.GT. jjm/2 +1 )  THEN
+       PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
+       STOP
+    ENDIF
+
+    IF( jfiltsu.LE.0) jfiltsu=1
+    IF( jfiltsu.GT.  jjm  +1 )  THEN
+       PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
+       STOP
+    ENDIF
+
+    IF( jfiltnv.LE.0) jfiltnv=1
+    IF( jfiltnv.GT. jjm/2    )  THEN
+       PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
+       STOP
+    ENDIF
+
+    IF( jfiltsv.LE.0) jfiltsv=1
+    IF( jfiltsv.GT.     jjm  )  THEN
+       PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
+       STOP
+    ENDIF
+
+    PRINT *,' jfiltnv jfiltsv jfiltnu jfiltsu ' , &
+         jfiltnv,jfiltsv,jfiltnu,jfiltsu
+
+    IF(first_call_inifilr) THEN
+       ALLOCATE(matriceun(iim,iim,jfiltnu))
+       ALLOCATE(matriceus(iim,iim,jfiltsu))
+       ALLOCATE(matricevn(iim,iim,jfiltnv))
+       ALLOCATE(matricevs(iim,iim,jfiltsv))
+       ALLOCATE( matrinvn(iim,iim,jfiltnu))
+       ALLOCATE( matrinvs(iim,iim,jfiltsu))
+       first_call_inifilr = .FALSE.
+    ENDIF
+
+    !                                 
+    !   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
+    !................................................................
+    !
+    !
+    DO j = 1,jjm
+       modfrstu( j ) = iim
+       modfrstv( j ) = iim
+    ENDDO
+    !
+    DO j = 2,jfiltnu
+       DO k = 2,modemax
+          cof = rlamda(k) * COS( rlatu(j) )
+          IF ( cof .LT. 1. ) GOTO 82
+       ENDDO
+       GOTO 84
+82     modfrstu( j ) = k
+       !
+       kf = modfrstu( j )
+       DO k = kf , modemax
+          cof = rlamda(k) * COS( rlatu(j) )
+          coefilu(k,j) = cof - 1.
+          coefilu2(k,j) = cof*cof - 1.
+       ENDDO
+84     CONTINUE
+    ENDDO
+    !                                 
+    !
+    DO j = 1,jfiltnv
+       !
+       DO k = 2,modemax
+          cof = rlamda(k) * COS( rlatv(j) )
+          IF ( cof .LT. 1. ) GOTO 87
+       ENDDO
+       GOTO 89
+87     modfrstv( j ) = k
+       !
+       kf = modfrstv( j )
+       DO k = kf , modemax
+          cof = rlamda(k) * COS( rlatv(j) )
+          coefilv(k,j) = cof - 1.
+          coefilv2(k,j) = cof*cof - 1.
+       ENDDO
+89     CONTINUE
+    ENDDO
+    !
+    DO j = jfiltsu,jjm
+       DO k = 2,modemax
+          cof = rlamda(k) * COS( rlatu(j) )
+          IF ( cof .LT. 1. ) GOTO 92
+       ENDDO
+       GOTO 94
+92     modfrstu( j ) = k
+       !
+       kf = modfrstu( j )
+       DO k = kf , modemax
+          cof = rlamda(k) * COS( rlatu(j) )
+          coefilu(k,j) = cof - 1.
+          coefilu2(k,j) = cof*cof - 1.
+       ENDDO
+94     CONTINUE
+    ENDDO
+    !                                 
+    DO j = jfiltsv,jjm
+       DO k = 2,modemax
+          cof = rlamda(k) * COS( rlatv(j) )
+          IF ( cof .LT. 1. ) GOTO 97
+       ENDDO
+       GOTO 99
+97     modfrstv( j ) = k
+       !
+       kf = modfrstv( j )
+       DO k = kf , modemax
+          cof = rlamda(k) * COS( rlatv(j) )
+          coefilv(k,j) = cof - 1.
+          coefilv2(k,j) = cof*cof - 1.
+       ENDDO
+99     CONTINUE
+    ENDDO
+    !
+
+    IF(jfiltnv.GE.jjm/2 .OR. jfiltnu.GE.jjm/2)THEN
+
+       IF(jfiltnv.EQ.jfiltsv)jfiltsv=1+jfiltnv
+       IF(jfiltnu.EQ.jfiltsu)jfiltsu=1+jfiltnu
+
+       PRINT *,'jfiltnv jfiltsv jfiltnu jfiltsu' , &
+            jfiltnv,jfiltsv,jfiltnu,jfiltsu
+    ENDIF
+
+    PRINT *,'   Modes premiers  v  '
+    PRINT 334,modfrstv
+    PRINT *,'   Modes premiers  u  '
+    PRINT 334,modfrstu
+
+    !  
+    !   ...................................................................
+    !
+    !   ... Calcul de la matrice filtre 'matriceu'  pour les champs situes
+    !                       sur la grille scalaire                 ........
+    !   ...................................................................
+    !
+    DO j = 2, jfiltnu
+
+       DO i=1,iim
+          coff = coefilu(i,j)
+          IF( i.LT.modfrstu(j) ) coff = 0.
+          DO k=1,iim
+             eignft(i,k) = eignfnv(k,i) * coff
+          ENDDO
+       ENDDO
+#ifdef CRAY
+       CALL MXM( eignfnv,iim,eignft,iim,matriceun(1,1,j),iim )
+#else
+#ifdef BLAS
+       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
+            eignfnv, iim, eignft, iim, 0.0, matriceun(1,1,j), iim)
+#else
+       DO k = 1, iim
+          DO i = 1, iim
+             matriceun(i,k,j) = 0.0
+             DO ii = 1, iim
+                matriceun(i,k,j) = matriceun(i,k,j) &
+                     + eignfnv(i,ii)*eignft(ii,k)
+             ENDDO
+          ENDDO
+       ENDDO
+#endif
+#endif
+
+    ENDDO
+
+    DO j = jfiltsu, jjm
+
+       DO i=1,iim
+          coff = coefilu(i,j)
+          IF( i.LT.modfrstu(j) ) coff = 0.
+          DO k=1,iim
+             eignft(i,k) = eignfnv(k,i) * coff
+          ENDDO
+       ENDDO
+#ifdef CRAY
+       CALL MXM(eignfnv,iim,eignft,iim,matriceus(1,1,j-jfiltsu+1),iim)
+#else
+#ifdef BLAS
+       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
+            eignfnv, iim, eignft, iim, 0.0, &
+            matriceus(1,1,j-jfiltsu+1), iim)
+#else
+       DO k = 1, iim
+          DO i = 1, iim
+             matriceus(i,k,j-jfiltsu+1) = 0.0
+             DO ii = 1, iim
+                matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) &
+                     + eignfnv(i,ii)*eignft(ii,k)
+             ENDDO
+          ENDDO
+       ENDDO
+#endif
+#endif
+
+    ENDDO
+
+    !   ...................................................................
+    !
+    !   ... Calcul de la matrice filtre 'matricev'  pour les champs situes
+    !                       sur la grille   de V ou de Z           ........
+    !   ...................................................................
+    !
+    DO j = 1, jfiltnv
+
+       DO i = 1, iim
+          coff = coefilv(i,j)
+          IF( i.LT.modfrstv(j) ) coff = 0.
+          DO k = 1, iim
+             eignft(i,k) = eignfnu(k,i) * coff
+          ENDDO
+       ENDDO
+#ifdef CRAY
+       CALL MXM( eignfnu,iim,eignft,iim,matricevn(1,1,j),iim )
+#else
+#ifdef BLAS
+       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
+            eignfnu, iim, eignft, iim, 0.0, matricevn(1,1,j), iim)
+#else
+       DO k = 1, iim
+          DO i = 1, iim
+             matricevn(i,k,j) = 0.0
+             DO ii = 1, iim
+                matricevn(i,k,j) = matricevn(i,k,j) &
+                     + eignfnu(i,ii)*eignft(ii,k)
+             ENDDO
+          ENDDO
+       ENDDO
+#endif
+#endif
+
+    ENDDO
+
+    DO j = jfiltsv, jjm
+
+       DO i = 1, iim
+          coff = coefilv(i,j)
+          IF( i.LT.modfrstv(j) ) coff = 0.
+          DO k = 1, iim
+             eignft(i,k) = eignfnu(k,i) * coff
+          ENDDO
+       ENDDO
+#ifdef CRAY
+       CALL MXM(eignfnu,iim,eignft,iim,matricevs(1,1,j-jfiltsv+1),iim)
+#else
+#ifdef BLAS
+       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
+            eignfnu, iim, eignft, iim, 0.0, & 
+            matricevs(1,1,j-jfiltsv+1), iim)
+#else
+       DO k = 1, iim
+          DO i = 1, iim
+             matricevs(i,k,j-jfiltsv+1) = 0.0
+             DO ii = 1, iim
+                matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) &
+                     + eignfnu(i,ii)*eignft(ii,k)
+             ENDDO
+          ENDDO
+       ENDDO
+#endif
+#endif
+
+    ENDDO
+
+    !   ...................................................................
+    !
+    !   ... Calcul de la matrice filtre 'matrinv'  pour les champs situes
+    !              sur la grille scalaire , pour le filtre inverse ........
+    !   ...................................................................
+    !
+    DO j = 2, jfiltnu
+
+       DO i = 1,iim
+          coff = coefilu(i,j)/ ( 1. + coefilu(i,j) )
+          IF( i.LT.modfrstu(j) ) coff = 0.
+          DO k=1,iim
+             eignft(i,k) = eignfnv(k,i) * coff
+          ENDDO
+       ENDDO
+#ifdef CRAY
+       CALL MXM( eignfnv,iim,eignft,iim,matrinvn(1,1,j),iim )
+#else
+#ifdef BLAS
+       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
+            eignfnv, iim, eignft, iim, 0.0, matrinvn(1,1,j), iim)
+#else
+       DO k = 1, iim
+          DO i = 1, iim
+             matrinvn(i,k,j) = 0.0
+             DO ii = 1, iim
+                matrinvn(i,k,j) = matrinvn(i,k,j) &
+                     + eignfnv(i,ii)*eignft(ii,k)
+             ENDDO
+          ENDDO
+       ENDDO
+#endif
+#endif
+
+    ENDDO
+
+    DO j = jfiltsu, jjm
+
+       DO i = 1,iim
+          coff = coefilu(i,j) / ( 1. + coefilu(i,j) )
+          IF( i.LT.modfrstu(j) ) coff = 0.
+          DO k=1,iim
+             eignft(i,k) = eignfnv(k,i) * coff
+          ENDDO
+       ENDDO
+#ifdef CRAY
+       CALL MXM(eignfnv,iim,eignft,iim,matrinvs(1,1,j-jfiltsu+1),iim)
+#else
+#ifdef BLAS
+       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
+            eignfnv, iim, eignft, iim, 0.0, matrinvs(1,1,j-jfiltsu+1), iim)
+#else
+       DO k = 1, iim
+          DO i = 1, iim
+             matrinvs(i,k,j-jfiltsu+1) = 0.0
+             DO ii = 1, iim
+                matrinvs(i,k,j-jfiltsu+1) = matrinvs(i,k,j-jfiltsu+1) &
+                     + eignfnv(i,ii)*eignft(ii,k)
+             ENDDO
+          ENDDO
+       ENDDO
+#endif
+#endif
+
+    ENDDO
+
+    !   ...................................................................
+
+    !
+334 FORMAT(1x,24i3)
+755 FORMAT(1x,6f10.3,i3)
+
+    RETURN
+  END SUBROUTINE inifilr
+
+END MODULE filtre
Index: LMDZ4/branches/LMDZ4-dev/libf/filtrez/parafilt.h
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/filtrez/parafilt.h	(revision 1074)
+++ LMDZ4/branches/LMDZ4-dev/libf/filtrez/parafilt.h	(revision 1086)
@@ -3,17 +3,2 @@
 !
         INTEGER nfilun, nfilus, nfilvn, nfilvs
-
-      PARAMETER (nfilun=24, nfilus=23, nfilvn=24, nfilvs=24)
-
-c
-c
-c      Ici , on a exagere  les nombres de lignes de latitudes a filtrer .
-c
-c      La premiere fois que  le Gcm  rentrera  dans le Filtre ,
-c
-c      il indiquera  les bonnes valeurs  de  nfilun , nflius, nfilvn  et 
-c
-c      nfilvs  a  mettre .  Il suffira alors de changer ces valeurs dans
-c
-c      Parameter  ci-dessus  et de relancer  le  run .  
-
