Index: trunk/LMDZ.MARS/libf/filtrez/filtreg_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/filtrez/inifilr.F	(revision 134)
+++ trunk/LMDZ.MARS/libf/filtrez/filtreg_mod.F90	(revision 1403)
@@ -1,16 +1,20 @@
+MODULE filtreg_mod
+
+CONTAINS
+
       SUBROUTINE inifilr
-c
-c    ... H. Upadhyaya, O.Sharma   ...
-c
+!
+!    ... H. Upadhyaya, O.Sharma   ...
+!
       IMPLICIT NONE
-c
-c     version 3 .....
-
-c     Correction  le 28/10/97    P. Le Van .
-c  -------------------------------------------------------------------
+!
+!     version 3 .....
+
+!     Correction  le 28/10/97    P. Le Van .
+!  -------------------------------------------------------------------
 #include "dimensions.h"
 #include "paramet.h"
 #include "parafilt.h"
-c  -------------------------------------------------------------------
+!  -------------------------------------------------------------------
 #include "comgeom.h"
 #include "coefils.h"
@@ -20,5 +24,5 @@
       REAL  dlonu(iim),dlatu(jjm)
       REAL  rlamda( iim ),  eignvl( iim )
-c
+!
 
       REAL    lamdamax,pi,cof
@@ -26,8 +30,8 @@
       REAL dymin,dxmin,colat0
       REAL eignft(iim,iim), coff
-      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 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)
 #ifdef CRAY
       INTEGER   ISMIN
@@ -37,23 +41,23 @@
 #endif
       EXTERNAL  inifgn
-c
-c ------------------------------------------------------------
-c   This routine computes the eigenfunctions of the laplacien
-c   on the stretched grid, and the filtering coefficients
-c      
-c  We designate:
-c   eignfn   eigenfunctions of the discrete laplacien
-c   eigenvl  eigenvalues
-c   jfiltn   indexof the last scalar line filtered in NH
-c   jfilts   index of the first line filtered in SH
-c   modfrst  index of the mode from where modes are filtered
-c   modemax  maximum number of modes ( im )
-c   coefil   filtering coefficients ( lamda_max*cos(rlat)/lamda )
-c   sdd      SQRT( dx )
-c      
-c     the modes are filtered from modfrst to modemax
-c      
-c-----------------------------------------------------------
-c
+!
+! ------------------------------------------------------------
+!   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. )
@@ -62,33 +66,33 @@
         dlonu(i) = xprimu( i )
        ENDDO
-c
+!
        CALL inifgn(eignvl)
-c
+!
         print *,' EIGNVL '
         PRINT 250,eignvl
 250     FORMAT( 1x,5e13.6)
-c
-c compute eigenvalues and eigenfunctions
-c
-c
-c.................................................................
-c
-c  compute the filtering coefficients for scalar lines and 
-c  meridional wind v-lines
-c
-c  we filter all those latitude lines where coefil < 1
-c  NO FILTERING AT POLES
-c
-c  colat0 is to be used  when alpha (stretching coefficient)
-c  is set equal to zero for the regular grid case 
-c
-c    .......   Calcul  de  colat0   .........
-c     .....  colat0 = minimum de ( 0.5, min dy/ min dx )   ...
-c
-c
+!
+! 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 45 j = 1,jjm
          dlatu( j ) = rlatu( j ) - rlatu( j+1 )
  45   CONTINUE
-c
+!
 #ifdef CRAY
       iymin   = ISMIN( jjm, dlatu, 1 )
@@ -106,30 +110,30 @@
        ENDDO
 #endif
-c
-c
+!
+!
       colat0  =  MIN( 0.5, dymin/dxmin )
-c
+!
       IF( .NOT.fxyhypb.AND.ysinus )  THEN
            colat0 = 0.6
-c         ...... a revoir  pour  ysinus !   .......
+!         ...... a revoir  pour  ysinus !   .......
            alphax = 0.
       ENDIF
-c
+!
       PRINT 50, colat0,alphax
   50  FORMAT(/15x,' Inifilr colat0 alphax ',2e16.7)
-c
+!
       IF(alphax.EQ.1. )  THEN
         PRINT *,' Inifilr  alphax doit etre  <  a 1.  Corriger '
          STOP
       ENDIF
-c
+!
       lamdamax = iim / ( pi * colat0 * ( 1. - alphax ) )
 
-cc                        ... Correction  le 28/10/97  ( P.Le Van ) ..
-c
+!c                        ... Correction  le 28/10/97  ( P.Le Van ) ..
+!
       DO 71 i = 2,iim
        rlamda( i ) = lamdamax/ SQRT( ABS( eignvl(i) ) )
  71   CONTINUE
-c
+!
 
       DO 72 j = 1,jjm
@@ -142,16 +146,16 @@
  72   CONTINUE
 
-c
-c    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....
-c    .........................................................
-c
+!
+!    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....
+!    .........................................................
+!
        modemax = iim
 
-cccc    imx = modemax - 4 * (modemax/iim)
+!ccc    imx = modemax - 4 * (modemax/iim)
 
        imx  = iim
-c
+!
        PRINT *,' TRUNCATION AT ',imx
-c
+!
       DO 75 j = 2, jjm/2+1
        cof = COS( rlatu(j) )/ colat0
@@ -162,9 +166,9 @@
        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
+          IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) &
+            jfiltsu= jjp1-j+1
         ENDIF
  75   CONTINUE
-c
+!
       DO 76 j = 1, jjm/2
        cof = COS( rlatv(j) )/ colat0
@@ -175,9 +179,9 @@
        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
+          IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) &
+             jfiltsv= jjm-j+1
         ENDIF
  76   CONTINUE
-c                                 
+!                                 
 
       IF( jfiltnu.LE.0 .OR. jfiltnu.GT. jjm/2 +1 )  THEN
@@ -201,17 +205,17 @@
       ENDIF
 
-       PRINT *,' jfiltnv jfiltsv jfiltnu jfiltsu ' ,
-     *           jfiltnv,jfiltsv,jfiltnu,jfiltsu
-
-c                                 
-c   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
-c................................................................
-c
-c
+       PRINT *,' jfiltnv jfiltsv jfiltnu jfiltsu ' , &
+                 jfiltnv,jfiltsv,jfiltnu,jfiltsu
+
+!                                 
+!   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
+!................................................................
+!
+!
       DO 77 j = 1,jjm
 	  modfrstu( j ) = iim
 	  modfrstv( j ) = iim
  77   CONTINUE
-c
+!
       DO 84 j = 2,jfiltnu
        DO 81 k = 2,modemax
@@ -221,5 +225,5 @@
       GOTO 84
  82   modfrstu( j ) = k
-c
+!
 	  kf = modfrstu( j )
 	   DO 83 k = kf , modemax
@@ -229,8 +233,8 @@
  83    CONTINUE
  84   CONTINUE
-c                                 
-c
+!                                 
+!
       DO 89 j = 1,jfiltnv
-c
+!
        DO 86 k = 2,modemax
 	    cof = rlamda(k) * COS( rlatv(j) )
@@ -239,5 +243,5 @@
       GOTO 89
  87   modfrstv( j ) = k
-c
+!
 	   kf = modfrstv( j )
 	   DO 88 k = kf , modemax
@@ -246,7 +250,7 @@
 	    coefilv2(k,j) = cof*cof - 1.
  88    CONTINUE
-c
+!
  89    CONTINUE
-c
+!
       DO 94 j = jfiltsu,jjm
        DO 91 k = 2,modemax
@@ -256,5 +260,5 @@
       GOTO 94
  92   modfrstu( j ) = k
-c
+!
         kf = modfrstu( j )
 	 DO 93 k = kf , modemax
@@ -264,5 +268,5 @@
  93      CONTINUE
  94    CONTINUE
-c                                 
+!                                 
       DO 99 j = jfiltsv,jjm
        DO 96 k = 2,modemax
@@ -272,5 +276,5 @@
       GOTO 99
  97   modfrstv( j ) = k
-c
+!
        kf = modfrstv( j )
 	   DO 98 k = kf , modemax
@@ -280,5 +284,5 @@
  98    CONTINUE
  99   CONTINUE
-c
+!
 
        IF(jfiltnv.GE.jjm/2 .OR. jfiltnu.GE.jjm/2)THEN
@@ -287,6 +291,6 @@
          IF(jfiltnu.EQ.jfiltsu)jfiltsu=1+jfiltnu
 
-          PRINT *,'jfiltnv jfiltsv jfiltnu jfiltsu' ,
-     *        jfiltnv,jfiltsv,jfiltnu,jfiltsu
+          PRINT *,'jfiltnv jfiltsv jfiltnu jfiltsu' , &
+              jfiltnv,jfiltsv,jfiltnu,jfiltsu
        ENDIF
 
@@ -298,87 +302,87 @@
      
       IF( nfilun.LT. jfiltnu )  THEN
-       PRINT *,' le parametre nfilun utilise pour la matrice ',
-     *   ' matriceun  est trop petit ! ' 
+       PRINT *,' le parametre nfilun utilise pour la matrice ', &
+         ' matriceun  est trop petit ! ' 
        PRINT *,'Le changer dans parafilt.h et le mettre a  ',jfiltnu
-        PRINT *,' Pour information, nfilun,nfilus,nfilvn,nfilvs '
-     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
-     *  ,jfiltnv,jjm-jfiltsv+1
+        PRINT *,' Pour information, nfilun,nfilus,nfilvn,nfilvs ' &
+       ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1 &
+        ,jfiltnv,jjm-jfiltsv+1
                STOP
       ENDIF
       IF( nfilun.GT. jfiltnu+ 2 )  THEN
-           PRINT *,' le parametre nfilun utilise pour la matrice ',
-     *' matriceun est trop grand ! Gachis de memoire ! ' 
+           PRINT *,' le parametre nfilun utilise pour la matrice ', &
+      ' matriceun est trop grand ! Gachis de memoire ! ' 
        PRINT *,'Le changer dans parafilt.h et le mettre a  ',jfiltnu
-        PRINT *,' Pour information, nfilun,nfilus,nfilvn,nfilvs '
-     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
-     *  ,jfiltnv,jjm-jfiltsv+1
-c              STOP
+        PRINT *,' Pour information, nfilun,nfilus,nfilvn,nfilvs ' &
+       ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1 &
+        ,jfiltnv,jjm-jfiltsv+1
+!              STOP
       ENDIF
       IF( nfilus.LT. jjm - jfiltsu +1 )  THEN
-            PRINT *,' le parametre nfilus utilise pour la matrice ',
-     *   ' matriceus  est trop petit !  '
-       PRINT *,' Le changer dans parafilt.h et le mettre a  ',
-     * jjm - jfiltsu + 1
-        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs '
-     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
-     *  ,jfiltnv,jjm-jfiltsv+1
+            PRINT *,' le parametre nfilus utilise pour la matrice ', &
+         ' matriceus  est trop petit !  '
+       PRINT *,' Le changer dans parafilt.h et le mettre a  ', &
+       jjm - jfiltsu + 1
+        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs ' &
+       ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1 &
+        ,jfiltnv,jjm-jfiltsv+1
                STOP
       ENDIF
       IF( nfilus.GT. jjm - jfiltsu + 3 )  THEN
-           PRINT *,' le parametre nfilus utilise pour la matrice ',
-     * ' matriceus  est trop grand ! ' 
-       PRINT *,' Le changer dans parafilt.h et le mettre a  ' ,
-     * jjm - jfiltsu + 1
-        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs '
-     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
-     *  ,jfiltnv,jjm-jfiltsv+1
-c              STOP
+           PRINT *,' le parametre nfilus utilise pour la matrice ', &
+       ' matriceus  est trop grand ! ' 
+       PRINT *,' Le changer dans parafilt.h et le mettre a  ' , &
+       jjm - jfiltsu + 1
+        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs ' &
+       ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1 &
+        ,jfiltnv,jjm-jfiltsv+1
+!              STOP
       ENDIF
       IF( nfilvn.LT. jfiltnv )  THEN
-            PRINT *,' le parametre nfilvn utilise pour la matrice ',
-     *   ' matricevn  est trop petit ! '  
+            PRINT *,' le parametre nfilvn utilise pour la matrice ', &
+         ' matricevn  est trop petit ! '  
        PRINT *,'Le changer dans parafilt.h et le mettre a  ',jfiltnv
-        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs '
-     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
-     *  ,jfiltnv,jjm-jfiltsv+1
+        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs ' &
+       ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1 &
+        ,jfiltnv,jjm-jfiltsv+1
                STOP
       ENDIF
       IF( nfilvn.GT. jfiltnv+ 2 )  THEN
-           PRINT *,' le parametre nfilvn utilise pour la matrice ',
-     *' matricevn est trop grand !  Gachis de memoire ! ' 
+           PRINT *,' le parametre nfilvn utilise pour la matrice ', &
+      ' matricevn est trop grand !  Gachis de memoire ! ' 
        PRINT *,'Le changer dans parafilt.h et le mettre a  ',jfiltnv
-        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs '
-     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
-     *  ,jfiltnv,jjm-jfiltsv+1
-c              STOP
+        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs ' &
+       ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1 &
+        ,jfiltnv,jjm-jfiltsv+1
+!              STOP
       ENDIF
       IF( nfilvs.LT. jjm - jfiltsv +1 )  THEN
-            PRINT *,' le parametre nfilvs utilise pour la matrice ',
-     *   ' matricevs  est trop petit !  Le changer dans parafilt.h '
-       PRINT *,' Le changer dans parafilt.h et le mettre a  '
-     * , jjm - jfiltsv + 1
-        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs '
-     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
-     *  ,jfiltnv,jjm-jfiltsv+1
+            PRINT *,' le parametre nfilvs utilise pour la matrice ', &
+         ' matricevs  est trop petit !  Le changer dans parafilt.h '
+       PRINT *,' Le changer dans parafilt.h et le mettre a  ' &
+       , jjm - jfiltsv + 1
+        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs ' &
+       ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1 &
+        ,jfiltnv,jjm-jfiltsv+1
                STOP
       ENDIF
       IF( nfilvs.GT. jjm - jfiltsv + 3 )  THEN
-           PRINT *,' le parametre nfilvs utilise pour la matrice ',
-     * ' matricevs  est trop grand ! Gachis de memoire ! '
-       PRINT *,' Le changer dans parafilt.h et le mettre a  '
-     *   ,  jjm - jfiltsv + 1
-        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs '
-     * ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1
-     *  ,jfiltnv,jjm-jfiltsv+1
-c              STOP
-      ENDIF
-
-c  
-c   ...................................................................
-c
-c   ... Calcul de la matrice filtre 'matriceu'  pour les champs situes
-c                       sur la grille scalaire                 ........
-c   ...................................................................
-c
+           PRINT *,' le parametre nfilvs utilise pour la matrice ', &
+       ' matricevs  est trop grand ! Gachis de memoire ! '
+       PRINT *,' Le changer dans parafilt.h et le mettre a  ' &
+         ,  jjm - jfiltsv + 1
+        PRINT *,' Pour information , nfilun,nfilus,nfilvn,nfilvs ' &
+       ,'doivent etre egaux successivement a  ',jfiltnu,jjm-jfiltsu+1 &
+        ,jfiltnv,jjm-jfiltsv+1
+!              STOP
+      ENDIF
+
+!  
+!   ...................................................................
+!
+!   ... Calcul de la matrice filtre 'matriceu'  pour les champs situes
+!                       sur la grille scalaire                 ........
+!   ...................................................................
+!
         DO j = 2, jfiltnu
 
@@ -394,6 +398,6 @@
 #else
 #ifdef BLAS
-         CALL SGEMM ('N', 'N', iim, iim, iim, 1.0,
-     $   eignfnv, iim, eignft, iim, 0.0, matriceun(1,1,j), iim)
+         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
@@ -401,6 +405,6 @@
             matriceun(i,k,j) = 0.0
             DO ii = 1, iim
-               matriceun(i,k,j) = matriceun(i,k,j)
-     .                          + eignfnv(i,ii)*eignft(ii,k)
+               matriceun(i,k,j) = matriceun(i,k,j) &
+                                + eignfnv(i,ii)*eignft(ii,k)
             ENDDO
          ENDDO
@@ -424,7 +428,7 @@
 #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)
+      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
@@ -432,6 +436,6 @@
             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)
+               matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) &
+                                          + eignfnv(i,ii)*eignft(ii,k)
             ENDDO
          ENDDO
@@ -442,10 +446,10 @@
         ENDDO
 
-c   ...................................................................
-c
-c   ... Calcul de la matrice filtre 'matricev'  pour les champs situes
-c                       sur la grille   de V ou de Z           ........
-c   ...................................................................
-c
+!   ...................................................................
+!
+!   ... Calcul de la matrice filtre 'matricev'  pour les champs situes
+!                       sur la grille   de V ou de Z           ........
+!   ...................................................................
+!
         DO j = 1, jfiltnv
 
@@ -461,6 +465,6 @@
 #else
 #ifdef BLAS
-      CALL SGEMM ('N', 'N', iim, iim, iim, 1.0,
-     $     eignfnu, iim, eignft, iim, 0.0, matricevn(1,1,j), iim)
+      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
@@ -468,6 +472,6 @@
             matricevn(i,k,j) = 0.0
             DO ii = 1, iim
-               matricevn(i,k,j) = matricevn(i,k,j)
-     .                          + eignfnu(i,ii)*eignft(ii,k)
+               matricevn(i,k,j) = matricevn(i,k,j) &
+                                + eignfnu(i,ii)*eignft(ii,k)
             ENDDO
          ENDDO
@@ -491,7 +495,7 @@
 #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)
+      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
@@ -499,6 +503,6 @@
             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)
+               matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) &
+                                + eignfnu(i,ii)*eignft(ii,k)
             ENDDO
          ENDDO
@@ -509,10 +513,10 @@
         ENDDO
 
-c   ...................................................................
-c
-c   ... Calcul de la matrice filtre 'matrinv'  pour les champs situes
-c              sur la grille scalaire , pour le filtre inverse ........
-c   ...................................................................
-c
+!   ...................................................................
+!
+!   ... Calcul de la matrice filtre 'matrinv'  pour les champs situes
+!              sur la grille scalaire , pour le filtre inverse ........
+!   ...................................................................
+!
         DO j = 2, jfiltnu
 
@@ -528,6 +532,6 @@
 #else
 #ifdef BLAS
-      CALL SGEMM ('N', 'N', iim, iim, iim, 1.0,
-     $     eignfnv, iim, eignft, iim, 0.0, matrinvn(1,1,j), iim)
+      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
@@ -535,6 +539,6 @@
             matrinvn(i,k,j) = 0.0
             DO ii = 1, iim
-               matrinvn(i,k,j) = matrinvn(i,k,j)
-     .                          + eignfnv(i,ii)*eignft(ii,k)
+               matrinvn(i,k,j) = matrinvn(i,k,j) &
+                                + eignfnv(i,ii)*eignft(ii,k)
             ENDDO
          ENDDO
@@ -558,6 +562,6 @@
 #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)
+      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
@@ -565,6 +569,6 @@
             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)
+               matrinvs(i,k,j-jfiltsu+1) = matrinvs(i,k,j-jfiltsu+1) &
+                                + eignfnv(i,ii)*eignft(ii,k)
             ENDDO
          ENDDO
@@ -575,10 +579,12 @@
         ENDDO
 
-c   ...................................................................
-
-c
+!   ...................................................................
+
+!
 334    FORMAT(1x,24i3)
 755    FORMAT(1x,6f10.3,i3)
 
        RETURN
-       END
+       END SUBROUTINE inifilr
+
+END MODULE  filtreg_mod
