--- trunk/libf/filtrez/inifilr.f90	2010/04/06 17:52:58	32
+++ trunk/Sources/filtrez/inifilr.f	2015/07/28 14:53:31	164
@@ -1,488 +1,262 @@
-SUBROUTINE inifilr
-
-  ! From filtrez/inifilr.F,v 1.1.1.1 2004/05/19 12:53:09
-  ! H. Upadhyaya, O.Sharma
-
-  !   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                    
-
-  USE dimens_m
-  USE paramet_m
-  USE logic
-  USE comgeom
-  USE serre
-  USE parafilt
-  USE coefils
+module inifilr_m
 
   IMPLICIT NONE
 
-  REAL dlonu(iim), dlatu(jjm)
-  REAL rlamda(iim), eignvl(iim)
+  ! North:
+
+  INTEGER jfiltnu, jfiltnv
+  ! index of the last scalar line filtered in northern hemisphere
 
-  REAL lamdamax, pi, cof
-  INTEGER i, j, modemax, imx, k, kf, ii
-  REAL dymin, dxmin, colat0
-  REAL eignft(iim,iim), coff
-  EXTERNAL inifgn
-
-  !-----------------------------------------------------------            
-
-
-  pi = 2.*asin(1.)
-
-  DO i = 1, iim
-     dlonu(i) = xprimu(i)
-  END DO
-
-  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)
-  END DO
-
-  dxmin = dlonu(1)
-  DO i = 2, iim
-     dxmin = min(dxmin,dlonu(i))
-  END DO
-  dymin = dlatu(1)
-  DO j = 2, jjm
-     dymin = min(dymin,dlatu(j))
-  END DO
-
-
-  colat0 = min(0.5,dymin/dxmin)
-
-  IF ( .NOT. fxyhypb .AND. ysinus) THEN
-     colat0 = 0.6
-     !         ...... a revoir  pour  ysinus !   .......                     
-     alphax = 0.
-  END IF
-
-  PRINT 50, colat0, alphax
-50 FORMAT (/15X,' Inifilr colat0 alphax ',2E16.7)
-
-  IF (alphax==1.) THEN
-     PRINT *, ' Inifilr  alphax doit etre  <  a 1.  Corriger '
-     STOP 1
-  END IF
-
-  lamdamax = iim/(pi*colat0*(1.-alphax))
-
-  DO  i = 2, iim
-     rlamda(i) = lamdamax/sqrt(abs(eignvl(i)))
-  END DO
-
-
-  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
-     end DO
-  END DO
-
-
-  !    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....          
-  !    .........................................................          
-
-  modemax = iim
-
-  !ccc    imx = modemax - 4 * (modemax/iim)                               
-
-  imx = iim
-
-  PRINT *, ' TRUNCATION AT ', imx
-
-  DO  j = 2, jjm/2 + 1
-     cof = cos(rlatu(j))/colat0
-     IF (cof<1.) THEN
-        IF (rlamda(imx)*cos(rlatu(j))<1.) jfiltnu = j
-     END IF
-
-     cof = cos(rlatu(jjp1-j+1))/colat0
-     IF (cof<1.) THEN
-        IF (rlamda(imx)*cos(rlatu(jjp1-j+1))<1.) jfiltsu = jjp1 - j + 1
-     END IF
-  END DO
-
-  DO  j = 1, jjm/2
-     cof = cos(rlatv(j))/colat0
-     IF (cof<1.) THEN
-        IF (rlamda(imx)*cos(rlatv(j))<1.) jfiltnv = j
-     END IF
-
-     cof = cos(rlatv(jjm-j+1))/colat0
-     IF (cof<1.) THEN
-        IF (rlamda(imx)*cos(rlatv(jjm-j+1))<1.) jfiltsv = jjm - j + 1
-     END IF
-  END DO
-
-
-  IF (jfiltnu<=0) jfiltnu = 1
-  IF (jfiltnu>jjm/2+1) THEN
-     PRINT *, ' jfiltnu en dehors des valeurs acceptables ', jfiltnu
-     STOP 1
-  END IF
-
-  IF (jfiltsu<=0) jfiltsu = 1
-  IF (jfiltsu>jjm+1) THEN
-     PRINT *, ' jfiltsu en dehors des valeurs acceptables ', jfiltsu
-     STOP 1
-  END IF
-
-  IF (jfiltnv<=0) jfiltnv = 1
-  IF (jfiltnv>jjm/2) THEN
-     PRINT *, ' jfiltnv en dehors des valeurs acceptables ', jfiltnv
-     STOP 1
-  END IF
-
-  IF (jfiltsv<=0) jfiltsv = 1
-  IF (jfiltsv>jjm) THEN
-     PRINT *, ' jfiltsv en dehors des valeurs acceptables ', jfiltsv
-     STOP 1
-  END IF
-
-  PRINT *, ' jfiltnv jfiltsv jfiltnu jfiltsu ', jfiltnv, jfiltsv, jfiltnu, &
-       jfiltsu
-
-
-  !   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....       
-  !................................................................       
-
-
-  DO  j = 1, jjm
-     modfrstu(j) = iim
-     modfrstv(j) = iim
-  END DO
-
-  DO  j = 2, jfiltnu
-     DO  k = 2, modemax
-        cof = rlamda(k)*cos(rlatu(j))
-        IF (cof<1.) GO TO 82
-     end DO
-     cycle
-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.
-     end DO
-  END DO
-
-
-  DO j = 1, jfiltnv
-     DO  k = 2, modemax
-        cof = rlamda(k)*cos(rlatv(j))
-        IF (cof<1.) GO TO 87
-     end DO
-     cycle
-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.
-     end DO
-  end DO
-
-  DO  j = jfiltsu, jjm
-     DO  k = 2, modemax
-        cof = rlamda(k)*cos(rlatu(j))
-        IF (cof<1.) GO TO 92
-     end DO
-     cycle
-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.
-     end DO
-  end DO
-
-  DO  j = jfiltsv, jjm
-     DO  k = 2, modemax
-        cof = rlamda(k)*cos(rlatv(j))
-        IF (cof<1.) GO TO 97
-     end DO
-     cycle
-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.
-     end DO
-  END DO
-
-
-  IF (jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2) THEN
-
-     IF (jfiltnv==jfiltsv) jfiltsv = 1 + jfiltnv
-     IF (jfiltnu==jfiltsu) jfiltsu = 1 + jfiltnu
-
-     PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', jfiltnv, jfiltsv, jfiltnu, &
-          jfiltsu
-  END IF
-
-  PRINT *, '   Modes premiers  v  '
-  PRINT 334, modfrstv
-  PRINT *, '   Modes premiers  u  '
-  PRINT 334, modfrstu
-
-
-  IF (nfilun<jfiltnu) THEN
-     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
-     STOP 1
-  END IF
-  IF (nfilun>jfiltnu+2) THEN
-     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
-  END IF
-  IF (nfilus<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
-     STOP 1
-  END IF
-  IF (nfilus>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
-  END IF
-  IF (nfilvn<jfiltnv) THEN
-     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
-     STOP 1
-  END IF
-  IF (nfilvn>jfiltnv+2) THEN
-     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
-  END IF
-  IF (nfilvs<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
-     STOP 1
-  END IF
-  IF (nfilvs>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
-  END IF
-
-  !   ... 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<modfrstu(j)) coff = 0.
-        DO k = 1, iim
-           eignft(i,k) = eignfnv(k,i)*coff
-        END DO
-     END DO
-     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)
-           END DO
-        END DO
-     END DO
-
-  END DO
-
-  DO j = jfiltsu, jjm
-
-     DO i = 1, iim
-        coff = coefilu(i,j)
-        IF (i<modfrstu(j)) coff = 0.
-        DO k = 1, iim
-           eignft(i,k) = eignfnv(k,i)*coff
-        END DO
-     END DO
-     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)
-           END DO
-        END DO
-     END DO
-
-  END DO
-
-  !   ................................................................... 
-
-  !   ... 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<modfrstv(j)) coff = 0.
-        DO k = 1, iim
-           eignft(i,k) = eignfnu(k,i)*coff
-        END DO
-     END DO
-     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)
-           END DO
-        END DO
-     END DO
-
-  END DO
-
-  DO j = jfiltsv, jjm
-
-     DO i = 1, iim
-        coff = coefilv(i,j)
-        IF (i<modfrstv(j)) coff = 0.
-        DO k = 1, iim
-           eignft(i,k) = eignfnu(k,i)*coff
-        END DO
-     END DO
-     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)
-           END DO
-        END DO
-     END DO
-
-  END DO
-
-  !   ................................................................... 
-
-  !   ... 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<modfrstu(j)) coff = 0.
-        DO k = 1, iim
-           eignft(i,k) = eignfnv(k,i)*coff
-        END DO
-     END DO
-     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)
-           END DO
-        END DO
-     END DO
-
-  END DO
-
-  DO j = jfiltsu, jjm
-
-     DO i = 1, iim
-        coff = coefilu(i,j)/(1.+coefilu(i,j))
-        IF (i<modfrstu(j)) coff = 0.
-        DO k = 1, iim
-           eignft(i,k) = eignfnv(k,i)*coff
-        END DO
-     END DO
-     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)
-           END DO
-        END DO
-     END DO
+  real, allocatable:: matriceun(:, :, :), matrinvn(:, :, :) 
+  ! (iim, iim, 2:jfiltnu)
 
-  END DO
+  real, allocatable:: matricevn(:, :, :) ! (iim, iim, jfiltnv)
+
+  ! South:
+
+  integer jfiltsu, jfiltsv
+  ! index of the first line filtered in southern hemisphere
+
+  real, allocatable:: matriceus(:, :, :), matrinvs(:, :, :)
+  ! (iim, iim, jfiltsu:jjm)
+
+  real, allocatable:: matricevs(:, :, :) ! (iim, iim, jfiltsv:jjm)
+
+contains
+
+  SUBROUTINE inifilr
+
+    ! From filtrez/inifilr.F, version 1.1.1.1, 2004/05/19 12:53:09
+    ! H. Upadhyaya, O. Sharma
+
+    ! This procedure computes the filtering coefficients for scalar
+    ! lines and meridional wind v lines. The modes are filtered from
+    ! modfrst to iim. 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.
+
+    USE dimens_m, ONLY : iim, jjm
+    USE dynetat0_m, ONLY : rlatu, rlatv, xprimu, grossismx
+    use inifgn_m, only: inifgn
+    use jumble, only: new_unit
+    use nr_util, only: pi
+
+    ! Local:
+
+    REAL dlatu(jjm)
+    REAL rlamda(2: iim)
+    real eignvl(iim) ! eigenvalues sorted in descending order (<= 0)
+    INTEGER i, j, unit
+    REAL colat0 ! > 0
+    REAL eignft(iim, iim)
+
+    real eignfnu(iim, iim), eignfnv(iim, iim)
+    ! eigenvectors of the discrete second derivative with respect to longitude
+
+    ! Filtering coefficients (lamda_max * cos(rlat) / lamda):
+    real, allocatable:: coefilnu(:, :) ! (iim, 2:jfiltnu)
+    real, allocatable:: coefilsu(:, :) ! (iim, jfiltsu:jjm)
+    real, allocatable:: coefilnv(:, :) ! (iim, jfiltnv)
+    real, allocatable:: coefilsv(:, :) ! (iim, jfiltsv:jjm)
+
+    ! Index of the mode from where modes are filtered:
+    integer, allocatable:: modfrstnu(:) ! (2:jfiltnu)
+    integer, allocatable:: modfrstsu(:) ! (jfiltsu:jjm)
+    integer, allocatable:: modfrstnv(:) ! (jfiltnv)
+    integer, allocatable:: modfrstsv(:) ! (jfiltsv:jjm)
+
+    !-----------------------------------------------------------
+
+    print *, "Call sequence information: inifilr"
+
+    CALL inifgn(eignvl, eignfnu, eignfnv)
+
+    ! Calcul de colat0
+    forall (j = 1:jjm) dlatu(j) = rlatu(j) - rlatu(j + 1)
+    colat0 = min(0.5, minval(dlatu) / minval(xprimu(:iim)))
+    PRINT *, 'colat0 = ', colat0
+
+    rlamda = iim / (pi * colat0 / grossismx) / sqrt(- eignvl(2: iim))
+
+    ! Determination de jfiltnu, jfiltsu, jfiltnv, jfiltsv
+
+    jfiltnu = (jjm + 1) / 2
+    do while (cos(rlatu(jfiltnu)) >= colat0 &
+         .or. rlamda(iim) * cos(rlatu(jfiltnu)) >= 1.)
+       jfiltnu = jfiltnu - 1
+    end do
+
+    jfiltsu = jjm / 2 + 2
+    do while (cos(rlatu(jfiltsu)) >= colat0 &
+         .or. rlamda(iim) * cos(rlatu(jfiltsu)) >= 1.)
+       jfiltsu = jfiltsu + 1
+    end do
+
+    jfiltnv = jjm / 2
+    do while ((cos(rlatv(jfiltnv)) >= colat0 &
+         .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) .and. jfiltnv >= 2)
+       jfiltnv = jfiltnv - 1
+    end do
+
+    if (cos(rlatv(jfiltnv)) >= colat0 &
+         .or. rlamda(iim) * cos(rlatv(jfiltnv)) >= 1.) then
+       ! {jfiltnv == 1}
+       PRINT *, 'Could not find jfiltnv.'
+       STOP 1
+    END IF
+
+    jfiltsv = (jjm + 1)/ 2 + 1
+    do while ((cos(rlatv(jfiltsv)) >= colat0 &
+         .or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) .and. jfiltsv <= jjm - 1)
+       jfiltsv = jfiltsv + 1
+    end do
+
+    IF (cos(rlatv(jfiltsv)) >= colat0 &
+         .or. rlamda(iim) * cos(rlatv(jfiltsv)) >= 1.) THEN
+       ! {jfiltsv == jjm}
+       PRINT *, 'Could not find jfiltsv.'
+       STOP 1
+    END IF
+
+    PRINT *, 'jfiltnu =', jfiltnu
+    PRINT *, 'jfiltsu =', jfiltsu
+    PRINT *, 'jfiltnv =', jfiltnv
+    PRINT *, 'jfiltsv =', jfiltsv
+
+    ! D\'etermination de coefil[ns][uv], modfrst[ns][uv]:
+
+    allocate(modfrstnu(2:jfiltnu), modfrstsu(jfiltsu:jjm))
+    allocate(modfrstnv(jfiltnv), modfrstsv(jfiltsv:jjm))
+    allocate(coefilnu(iim, 2:jfiltnu), coefilsu(iim, jfiltsu:jjm))
+    allocate(coefilnv(iim, jfiltnv), coefilsv(iim, jfiltsv:jjm))
+
+    coefilnu = 0.
+    coefilnv = 0.
+    coefilsu = 0.
+    coefilsv = 0.
+
+    DO j = 2, jfiltnu
+       modfrstnu(j) = 2
+       do while (rlamda(modfrstnu(j)) * cos(rlatu(j)) >= 1. &
+            .and. modfrstnu(j) <= iim - 1)
+          modfrstnu(j) = modfrstnu(j) + 1
+       end do
+
+       if (rlamda(modfrstnu(j)) * cos(rlatu(j)) < 1.) then
+          DO i = modfrstnu(j), iim
+             coefilnu(i, j) = rlamda(i) * cos(rlatu(j)) - 1.
+          end DO
+       end if
+    END DO
+
+    DO j = 1, jfiltnv
+       modfrstnv(j) = 2
+       do while (rlamda(modfrstnv(j)) * cos(rlatv(j)) >= 1. &
+            .and. modfrstnv(j) <= iim - 1)
+          modfrstnv(j) = modfrstnv(j) + 1
+       end do
+
+       if (rlamda(modfrstnv(j)) * cos(rlatv(j)) < 1.) then
+          DO i = modfrstnv(j), iim
+             coefilnv(i, j) = rlamda(i) * cos(rlatv(j)) - 1.
+          end DO
+       end if
+    end DO
+
+    DO j = jfiltsu, jjm
+       modfrstsu(j) = 2
+       do while (rlamda(modfrstsu(j)) * cos(rlatu(j)) >= 1. &
+            .and. modfrstsu(j) <= iim - 1)
+          modfrstsu(j) = modfrstsu(j) + 1
+       end do
+
+       if (rlamda(modfrstsu(j)) * cos(rlatu(j)) < 1.) then
+          DO i = modfrstsu(j), iim
+             coefilsu(i, j) = rlamda(i) * cos(rlatu(j)) - 1.
+          end DO
+       end if
+    end DO
+
+    DO j = jfiltsv, jjm
+       modfrstsv(j) = 2
+       do while (rlamda(modfrstsv(j)) * cos(rlatv(j)) >= 1. &
+            .and. modfrstsv(j) <= iim - 1)
+          modfrstsv(j) = modfrstsv(j) + 1
+       end do
+
+       if (rlamda(modfrstsv(j)) * cos(rlatv(j)) < 1.) then
+          DO i = modfrstsv(j), iim
+             coefilsv(i, j) = rlamda(i) * cos(rlatv(j)) - 1.
+          end DO
+       end if
+    END DO
+
+    call new_unit(unit)
+    open(unit, file = "inifilr_out.txt", status = "replace", action = "write") 
+    write(unit, fmt = *) '"EIGNVL"', eignvl
+    write(unit, fmt = *) '"modfrstnu"', modfrstnu
+    write(unit, fmt = *) '"modfrstsu"', modfrstsu
+    write(unit, fmt = *) '"modfrstnv"', modfrstnv
+    write(unit, fmt = *) '"modfrstsv"', modfrstsv
+    close(unit)
+
+    allocate(matriceun(iim, iim, 2:jfiltnu), matrinvn(iim, iim, 2:jfiltnu))
+    allocate(matricevn(iim, iim, jfiltnv))
+    allocate(matricevs(iim, iim, jfiltsv:jjm))
+    allocate(matriceus(iim, iim, jfiltsu:jjm), matrinvs(iim, iim, jfiltsu:jjm))
+
+    ! Calcul de la matrice filtre 'matriceu' pour les champs situes
+    ! sur la grille scalaire
+
+    DO j = 2, jfiltnu
+       eignft(:modfrstnu(j) - 1, :) = 0.
+       forall (i = modfrstnu(j):iim) eignft(i, :) = eignfnv(:, i) &
+            * coefilnu(i, j)
+       matriceun(:, :, j) = matmul(eignfnv, eignft)
+    END DO
+
+    DO j = jfiltsu, jjm
+       eignft(:modfrstsu(j) - 1, :) = 0.
+       forall (i = modfrstsu(j):iim) eignft(i, :) = eignfnv(:, i) &
+            * coefilsu(i, j)
+       matriceus(:, :, j) = matmul(eignfnv, eignft)
+    END DO
+
+    ! Calcul de la matrice filtre 'matricev' pour les champs situes
+    ! sur la grille de V ou de Z
+
+    DO j = 1, jfiltnv
+       eignft(:modfrstnv(j) - 1, :) = 0.
+       forall (i = modfrstnv(j): iim) eignft(i, :) = eignfnu(:, i) &
+            * coefilnv(i, j)
+       matricevn(:, :, j) = matmul(eignfnu, eignft)
+    END DO
+
+    DO j = jfiltsv, jjm
+       eignft(:modfrstsv(j) - 1, :) = 0.
+       forall (i = modfrstsv(j):iim) eignft(i, :) = eignfnu(:, i) &
+            * coefilsv(i, j)
+       matricevs(:, :, j) = matmul(eignfnu, eignft)
+    END DO
+
+    ! Calcul de la matrice filtre 'matrinv' pour les champs situes
+    ! sur la grille scalaire , pour le filtre inverse
+
+    DO j = 2, jfiltnu
+       eignft(:modfrstnu(j) - 1, :) = 0.
+       forall (i = modfrstnu(j):iim) eignft(i, :) = eignfnv(:, i) &
+            * coefilnu(i, j) / (1. + coefilnu(i, j))
+       matrinvn(:, :, j) = matmul(eignfnv, eignft)
+    END DO
+
+    DO j = jfiltsu, jjm
+       eignft(:modfrstsu(j) - 1, :) = 0.
+       forall (i = modfrstsu(j):iim) eignft(i, :) = eignfnv(:, i) &
+            * coefilsu(i, j) / (1. + coefilsu(i, j))
+       matrinvs(:, :, j) = matmul(eignfnv, eignft)
+    END DO
 
-334 FORMAT (1X,24I3)
-755 FORMAT (1X,6F10.3,I3)
+  END SUBROUTINE inifilr
 
-END SUBROUTINE inifilr
+end module inifilr_m