Index: /LMDZ4/branches/LMDZ4-dev/libf/filtrez/filtreg_mod.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/filtrez/filtreg_mod.F90	(revision 1109)
+++ /LMDZ4/branches/LMDZ4-dev/libf/filtrez/filtreg_mod.F90	(revision 1109)
@@ -0,0 +1,530 @@
+MODULE filtreg_mod
+
+  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 filtreg_mod
Index: DZ4/branches/LMDZ4-dev/libf/filtrez/module_filtreg.F90
===================================================================
--- /LMDZ4/branches/LMDZ4-dev/libf/filtrez/module_filtreg.F90	(revision 1108)
+++ 	(revision )
@@ -1,530 +1,0 @@
-MODULE filtreg_mod
-
-  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 filtreg_mod